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Sampling From Structured Populations: Some 
Issues and Answers 


By V. N. NAIR and T. E. DALENIUS* 
(Manuscript received February 25, 1981) 


This paper reviews some sampling issues that are common to many 
Bell System surveys. We discuss various aspects of two-stage sam-. 
pling designs, and emphasize sampling from populations with mul- 
tiple characteristics. The hierarchical structure of the population in 
many surveys makes the use of multistage sampling techniques at- 
tractive. In populations with multiple characteristics, often not every 
characteristic is common to every unit. We consider some special 
designs for sampling from such populations. Finally, we discuss some 
issues in network sampling. Two recent Bell System surveys are used 
to illustrate most of the ideas discussed. One of the surveys deals with 
the estimation of traffic characteristics for various classes of service, 
while the other one is a survey of baseband transmission impairments. 


I. INTRODUCTION 


Sample surveys have played an increasingly important role in the 
Bell System in recent years as a means of providing an objective basis 
for decision making. To an extent, this has been due to the growing 
awareness among users of the survey results that, in most surveys, 
sampling is not the only source of error and often not the primary 
source. Even if a presumably complete census were taken instead of a 
sample, serious errors might exist in the results arising from various 
causes such as measurement or response errors. 


* Brown University. 
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The growth in numbers, in recent years, has also been accompanied 
by a widening of the range (both in type and complexity) of the 
surveys. For many of these surveys, a simple and readily available 
sampling design can easily be adapted to the needs of the prevailing 
situation. More often, however, the problem at hand is sufficiently 
complex and nonstandard so that various parts of existing sampling 
theory have to be modified and pieced together to arrive at a reason- 
able solution. 

Nevertheless, some sampling issues are common to a number of Bell 
System surveys. Most of these surveys involve sampling from popula- 
tions that are highly structured, and any cost-efficient sampling design 
must take this structure into account. In this paper, we review some 
sampling issues that arose in two surveys currently under implemen- 
tation. Both surveys possess some common features as well as features 
unique to themselves. Since these features are common to a large 
number of other surveys, an exposition of both the theoretical and 
practical considerations involved may prove beneficial to other survey 
practitioners. Let us first consider the two examples. 


Example 1. Cost of service traffic usage studies (COSTUS) 


The various Bell operating telephone companies (OTCs) carry out 
these surveys periodically to obtain an objective basis for distributing 
the traffic-sensitive costs for a jurisdiction, typically a state within an 
OTC, among its various classes of telephone service. Measurements of 
three traffic characteristics (busy-hour ccs, busy-hour peg count and 
14-day peg count) from the sampled telephone lines are used to 
calculate the relative magnitudes of the traffic characteristics for each 
class of service. [CCS is a traditional unit for measuring the usage of 
channels (it stands for hundred call seconds per hour). Peg count is 
the number of calls actually handled.] These values are then used as 
inputs to the “embedded direct costs” analysis, which allocates most 
traffic-sensitive investments and expenses among the various classes 
of service. 

The elementary units in this study are telephone lines corresponding 
to the various classes of service. These units, however, are clustered 
into central offices. In fact, each central office has a number of clusters 
associated with it, one cluster for each class of service. A reasonably 
cost-efficient design should take this hierarchical clustering into ac- 
count, since the major portion of the costs in observing a line arises 
from visiting the central office and setting up the measuring equip- 
ment. Thus, a two-stage sampling design with central offices serving 
as primary sampling units (psUs) and telephone lines serving as sec- 
ondary sampling units (ssus) seems attractive. This is even more so 
since the central offices provide service in a number of classes of 
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service so that from each sampled central office, we can further 
subsample telephone lines from all the available classes of service. 
Hence, cosTUS are examples of the use of a two-stage sampling 
design for a population with multiple characteristics. The different 
characteristics here correspond to the different classes of service. The 
parameters of the sampling design in cosTUS are determined so that 
the busy-hour ccs parameter for each class is estimated with a pre- 
scribed accuracy. One additional complication in these studies is the 
fact that not all central offices provide service in every available class. 
In some jurisdictions, there are some classes of service (such as coin) 
that are provided in only a few offices. The sampling literature refers 
to this as the problem of “partial variate pattern” (pvp). The presence 
of PvP causes difficulties in selecting an appropriate sample of central 
offices for the estimation of the parameters of all the classes of service. 


Example 2. Survey of baseband transmission impairments 


The aim of this survey, currently under development at Bell Labo- 
ratories, is to measure baseband transmission impairments for various 
trunk facility types. From each sampled trunk, estimates of various 
impairment characteristics, such as signal to C-notched noise ratio (s/ 
n) and second- and third-order harmonic distortion (R2 and R3) are to 
be obtained. Although the near (transmitting) and far (receiving) end- 
drop equipment, in addition to the carrier system, determines the 
trunk type, it is known from past experience that the contribution 
from the carrier system is the dominant factor. Thus, we do not 
consider the influence of the end-drop equipment in this study. Six 
different measurement characteristics are to be measured from each 
sampled trunk and the parameters of seven different trunk types are 
to be estimated. 

The elementary unit in this survey is the trunk. While the trunks 
are again clustered into central offices, this clustering is not unique 
since one trunk is common to a pair (transmitting and receiving) of 
central offices. In fact, the structure of the population here resembles 
a graph (network) with the central offices as nodes and trunks as edges 
(arcs). This survey is an example of network (graph) sampling (see 
Ref. 1, for example). In this survey, if we sample a particular trunk, we 
have to visit the pair of end offices connected to the trunk to set up 
the measuring equipment. This implies that it is cheaper to sample 
additional trunks connected to those two end offices. Hence, taking 
the structure of the population into account results in considerable 
cost savings. . 

One possible approach to this problem is to use multistage sampling 
to select pairs of offices and trunks connected to those offices. Since 
we are interested in different trunk types, this study also involves 
multiple characteristics. 
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Both the above examples involve using multistage sampling to study 
populations with multiple characteristics. Multistage sampling is not 
an uncommon phenomenon in Bell System surveys where the natural 
administrative and geographic clustering of units makes it very cost 
efficient. In Sections II and III we review various issues that confront 
a survey statistician in developing a two-stage sampling design for 
studying multiple characteristics. Some of the issues discussed in 
Section II are also common to other sampling designs. Section III 
deals primarily with determining the parameters of the sample design. 
In Section IV, we consider some sampling designs for populations with 
Pvp. Section V is a brief review of issues in network sampling. We 
conclude the paper with a summary in Section VI. Throughout the 
paper we try to balance theoretical considerations with practical 
guidelines gained from our own experience. One of the two examples 
is used, wherever possible, to illustrate the ideas discussed. 


ll. TWO-STAGE SAMPLING: SOME PRELIMINARIES 


This section deals with some preliminary considerations in devel- 
oping a two-stage sampling design. Some of the discussion deals with 
issues that are common to sample surveys in general. We begin with 
a discussion of the rationale for using two-stage or multistage sampling 
designs. After an introduction to some notation, we examine how 
prescribed accuracy requirements are implemented in a sample survey 
and discuss the use of prior information. Section 2.6 examines the use 
of varying probability sampling schemes. Section 2.7 discusses ratio 
estimators with specific emphasis on two-stage sampling situations. 


2.1 Why two-stage sampling? 


The individuals whose characteristics are to be measured in a study 
are called elementary units. Observational access to the elementary 
units, in many cases, is provided by multistage sampling. Let the 
elementary units be grouped into a number of suitable clusters. In two- 
stage sampling, the clusters are used as PSUs and a sample of Psus is 
selected in the first stage. The Psus selected are divided into a number 
of ssUs and a sample of ssus is selected from each Psu selected in the 
first stage. (The elementary units themselves can serve as ssus.) All 
elementary units in the selected ssus are observed with respect to the 
variables of interest. 

There are various reasons why multistage sampling is attractive. For 
instance, in many studies, a complete list (“frame”) of elementary 
units is not available and it may be prohibitively expensive to create 
such a list. If it is relatively cheap to construct a list of clusters, the 
clusters can be used as PSUs in a two-stage sampling scheme. Then, 
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only a list of the elementary units in the sampled clusters needs to be 
constructed. This results in considerable cost savings. 

Often, the population of elementary units in a survey is dispersed 
over a large geographical area. If we have to visit each sampled unit to 
collect measurements, sampling from the list of elementary units can 
lead to high costs per elementary unit. A more cost-efficient scheme 
may be obtained by grouping the elementary units into geographically 
compact clusters and using multistage sampling with the clusters as 
PSUS. 

Typically, the cost reduction in multistage sampling is accompanied 
by an increase in the variance of the estimate over the variance of an 
estimate from a simple random sampling (SRS) of the same number of 
elementary units. However, the “accuracy” per unit cost may be 
higher. If we have some control over the formation of the clusters, we 
can actually reduce the variance (relative to SRS) by grouping the units 
so that there is more variation within clusters than between clusters. 
In most Bell System surveys, however, the clusters are predetermined. 


2.2 Notation 


We use the following notation throughout the remainder of this 
paper: 


M = number of psus in the universe, 


m = number of psus sampled, 

N; = number of ssUsin PSUI, 1=1,---,M, 

n; = number of ssus selected from the ith sampled psu, 
L=1,--+,m, 

II; = probability of selecting the ith Psu in a sample of size 


Mm, ye) Il; =m, 
Y,; = characteristic to be measured, 7=1,---,Ni, t=1,---, 


M, 
yi = value corresponding to a sample unit, j=1,--- , ni, 
1 = I, eee, m, 
N; M = 
Y; = > Yi; Y= > Yi, Y; = Y;/Ni, 
= ie M 1 Ni = 
Y = Y/N, N=% Ni, S?=— ¥ (Y;—- Yi)’, 
i=1 N; j=1 


Des Ves » Jijs y=) % Vi = yi/Ni, 
= 
y = y/n, n= >; Ni. 


We consider only equal probability sampling schemes in stage 
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two in this paper. The parameter of interest is the overall total Y = 
y, Ni¥;. ¥ denotes an arbitrary estimator of Y. The same consider- 
ations can be used for estimating the average Y if we rewrite 


_ M = 
=) W.Y,, W; = Ni/N. 


i=1 


2.3 Accuracy requirements 


The sampling design in a carefully planned survey is determined so 
that either (z) the total cost of the survey is minimized subject to a 
prescribed requirement on the accuracy of the estimators or (ii) the 
accuracy of the estimators is maximized subject to a constraint on the 
cost. Since both approaches involve essentially the same considerations 
(see Section III), let us consider in some detail just the problem of 
minimizing cost subject to accuracy requirements. 

A sampling design, where the units are randomly selected according 
to given probabilities of selection, permits us to make quantitative 
statements about the error involved in the estimators. This in turn 
allows us to determine the sample sizes so that the prescribed accuracy 
requirements are met. These requirements are typically stated in terms 
of the error e = Y — Y or some function of the error, f(e), such as 
relative error, and can be expressed as 


Pr{| f(e)|<d}=1—@ (1) 


for some constants a and 6. In costus, for instance, the sample sizes 
are determined so that the absolute values of the relative error is less 
than or equal to 0.1 with probability at least 0.9, i.e., a = 6 = 0.1. To 
implement the accuracy condition (1), large-sample theory is usually 
used to claim that Y is approximately normally distributed. (It is 
beyond the scope of this paper to discuss the adequacy of this normal 
approximation. The interested reader is referred to Refs. 2 to 5.) 
Equation (1) is equivalent to an expression of an upper | bound on the 
variance [or mean-square error (mse) if ¥ is biased] of Y. 

When estimating several parameters, as in populations with multiple 
characteristics, we may require that several accuracy criteria be sat- 
isfied simultaneously. By using normal approximations, we can state 
this problem, in general, as minimizing’ the total cost of the survey 
subject to a constraint on the variances (or mse’s) of the form 


Av =}, (2) 


where v = (Uj, «++, Up)” is the vector of variances (mse’s) of the p 
estimators, A is a k X p matrix that specifies the k specific linear 
combinations of the variances that have to meet the accuracy condi- 
tions, and y = (y, «++ , yx)” represents the bounds on the accuracies. 
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For example, if k = p and A is the identity matrix, then all the p 
parameters need to be estimated with prescribed accuracy. If k = 1, 
then only one particular linear combination of the variances is needed 
to satisfy an accuracy criterion. 


2.4 Variance components 


Since the accuracy specifications can be stated in terms of the 
variances of the individual estimators, we need to examine the com- 
ponents of the variance of the estimator in a two-stage sampling 
scheme. This will aid us later (Section III) in determining the relative 
contributions to the variance from stages one and two and the tradeoffs 
in increasing the sample size in stage one versus that in stage two. If 
we restrict our attention to linear estimators of the form Y., = Viaiyi 
for estimating Y, we see that a; must equal N;,/TI; for the estimator to 
be unbiased. With this choice of a, Y,, is the well-known Horvitz- 
Thompson (H-T) estimator.® A discussion of some of the properties of 
this estimator can be found in Ref. 7. Let us restrict our attention to 
the H-T estimator and examine its variance. 

If we select m Psus with replacement (WR) in stage one with inclusion 
probabilities IT;, we have a multinomial sample of size m with success 
probabilities Z; = I1;/m. If the second-stage units are chosen without 
replacement, the variance of 


can be written as the sum of two components:”® 
(1) the within-Psu variation W is 
1X NS 
> Zn (1 — fai), (3) 


W=— 
M j=1 Li 


and 
(ti) the between-Psu variation B is 


1 M 
B=—Y ZiA(Yi/Z: — Y)’. 
M j=1 


Here, 


the within cluster variance and 1 — fo; = (Ni — nj)/Ni, the finite 
population correction. 

If the sampling is done without replacement (wor) in stage one with 
varying selection probabilities, the within-psu variation remains the 
same. The between-psu variation, however, depends on second-order 
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inclusion probabilities which are extremely hard to calculate.’’° Har- 
tley and Rao provide some approximations.” One possible approxi- 
mation is, of course, the use of eq. (3), valid for the WR scheme, in the 
WOR Situation. If the sampling fraction m/M is large (say >0.25), this 
approximation may be unreasonable. When the sampling is done wor 
with equal selection probabilities in stage one, 1.e., SRSwoR, the B 
component is given by 


_M(1-f) S _ v2 
ie eo 


where Y = a yw, Y; and f= m/M. 


For a discussion of variance estimation in two-stage sampling, see 
Refs. 7, 8, or 9, for example. Some approximate but “quick and easy” 
methods of variance estimation are discussed in Refs. 11 and 12. If the 
variance estimator is intended only to provide a rough guide as to the 
accuracy of the estimator, an approximate, but quick and easy, method 
is adequate. If the accuracy of the estimator is of great importance 
and must be demonstrated through the variance estimator, we have to 
use a “good” variance estimator, such as one with small mse. 


2.5 Prior information 


We need prior information on the variance of the various estimators 
and on the sampling costs to determine the sample sizes in a survey. 
It is rare that we have very good prior information, particularly 
concerning the variance of the estimators. Preliminary estimates can 
be obtained from prior surveys or pilot studies. One practice commonly 
found in the Bell System is the use of data from the entire Bell System 
to develop preliminary estimates for specific jurisdictions. 

To implement the accuracy conditions exactly in a two-stage sam- 
pling scheme, we need to know each one of the components of W and 
B in eq. (3) exactly. Since this is rather unlikely, we usually just use 
two numbers, one for W and one for B, instead of the individual values 
for each psu. These numbers can be interpreted as either the average 
or the maximum over all Psus. 

When the quality of the prior information is poor (as a consequence 
of one or more of the above reasons), little can be gained in developing 
a complex design that may (or may not) be “optimum” for the problem 
at hand. A simpler design which is less sensitive to the preliminary 
estimates of the design parameters is more desirable. Also, when the 
preliminary variance estimators are unreliable, an estimate of the 
accuracy achieved should always be calculated after the fact from the 
sample to compare with the prescribed accuracy. 
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2.6 Varying probability sampling 


The sample selection schemes in stages one and two can be based 
on equal or varying probability sampling techniques. For simplicity, 
we consider varying probability sampling only in stage one. The 
considerations here also carry over to other stages. Let us examine 
how the selection probabilities {II;} should be determined so that the 
variance of the H-T estimator 5”, (N;/II,)y;, for estimating Y = 
yt, Y,, is minimized. 

In the simpler situation of one-stage cluster sampling, i.e., n; = Ni, 
if we take IT; proportional to Y;, the variance of the H-T estimator is 
zero.’ Hence, if there exists an auxiliary variable X; which is approxi- 
mately proportional to Y;, we can use this auxiliary information to 
select the IT,’s. In some two-stage sampling situations, we can use the 
measures of size of the psu, {Ni}, to obtain “optimal” selection prob- 
abilities. To see this, note that the parameter Y can be written as 
yi, N:Yi, where Y; is the psu mean, and often the Y;’s are roughly of 
the same order of magnitude. In this case, the Y = NY; will be roughly 
proportional to the N; so that we can take the II; proportional to the 
N;. This is known as probability proportional to size (PPS) sampling. 
(In cosTus, for example, a priori, we expect the average busy-hour ccs 
per main station to be about the same across central offices.) When 
sampling from populations with multiple characteristics, there are 
multiple measures of size, one associated with each characteristic. The 
optimal selection probabilities are some function of these size mea- 
sures, depending on the particular accuracy criteria of interest. In 
addition, there are also cases in which the exact size measures are 
unknown and we have to use estimated measures. . 

To develop a cost-efficient design, we need to minimize variance per 
unit cost rather than the actual variance. The optimal selection prob- 
abilities must therefore take the cost structure into account. In COSTUS, 
where the psus are central offices, the sampling costs depend on the 
type of switching equipment in the office. For example, it is consider- 
ably more expensive to visit and set up the measuring equipment in an 
electronic switching system (Ess) office than in a non-Ess office. If we 
use formal optimality calculations, we find that with other factors held 
constant, the optimal selection probability for each Psu is inversely 
proportional to the square root of the cost of sampling that psu.® 

One or more of the above considerations may indicate that even if 
the PsUs vary greatly in size, the optimal selection probabilities are 
not too unequal. In such a case, we may be better off using SRS, Le., 
equal selection probabilities, since (i) the selection scheme is simpler 
and (iz) exact variance formulas are available if, in addition, we are 
sampling wor. In some situations, we can actually calculate the gain 
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from using varying selection probability schemes.’ If the gain is not 
substantial in these situations, the use of SRS seems preferable. 

Also, even if we use SRS when the psus vary greatly in size, we can 
use ratio estimators, which take into account this variation, to estimate 
the parameters. This is discussed in Section 2.7. 

Finally, we briefly discuss a simple scheme for selecting PsuUs with 
unequal probabilities. Many schemes for unequal probability selection 
exist,” and, in fact, several procedures may lead to the same inclusion 
probabilities {II;}. The scheme we consider here is for sampling woR 
and is known as PPS systematic sampling. Let {7;} denote the cumu- 
lative totals of the ape selection probabilities {TI}, 


5 II; =m, T; = y II. 

i=] j=1 
To select m psus, first select a random number u € [0, 1] and then 
select the m psus for which 


Ti1<ut+js Ti, J=0,1,---,m-—1. 


Hartley and Rao consider this procedure with a random arrangement 
of the psus and develop approximate variance expressions for the 
estimator.’® 


2.7 Use of ratio estimation 


So far we have considered only unbiased estimators of the total Y. 
In some situations we can exploit information available for some 
auxiliary variable and use a biased estimator, such as the ratio esti- 
mator, which has smaller mse than the unbiased estimators. To see 
this, let {X;} be the known auxiliary variable and let X denote the total 
corresponding to this variable and X denote the estimator of X based 
on the sample. Since we know the error X — X, we know how this 
sample performs in estimating X. Hence, if {X;} and {Y,} are highly 
correlated, it is intuitively clear that we can improve our original 
estimator Y by exploiting our knowledge of how well the sample 
estimates X. 

The ratio estimator itself is a special case of the general difference 
estimator Y, = Y + a(X — X) and i is obtained by taking a = —Y¥/X. 
This results in the estimator Y = (Y/X )X. There are other ways of 
exploiting the information about X — X. For instance, a can be a 
prespecified constant. (If a = 0, we get the original estimator Y based 
on the Y measurements alone.) We can also take a to be the regression 
coefficient B obtained by regressing the Y;’s on the X’s. 

For the ratio estimator Y, the mse of Y can be approximated up to 
a first-order term by 


1244 THE BELL SYSTEM TECHNICAL JOURNAL, SEPTEMBER 1981 


V(Y) — 2R Cov(Y, X) + R?V(X), 


where R = Y/X and Y and X are unbiased estimators of Y and X.° 
Thus, Y will be more efficient than the unbiased estimator Y if 
2R Cov(Y, X) > R2V(X). This is likely to be true in practice if the X/’s 
are appropriately chosen. 

In Section 2.6, we saw that in some two-stage sampling situations, 
the Y;,’s are likely to be correlated with the size measures {N;}. If pps 
sampling is not used (for one or more of the reasons we considered 
earlier), we can take the N;’s to be the auxiliary variables and use the 
resulting ratio estimator Y = RN, where 


a Ni 5 a 
R= 47? ue — 
(If we use PPS sampling, the ratio estimator with the size measure as 
the auxiliary variable is the same as the unbiased estimator.) Our 


experience with data from several jurisdictions for cosTUS showed a 
considerable gain from the use of this ratio estimator. 


Ii!! DETERMINING THE DESIGN PARAMETERS 
3.1 Cost considerations 


The ultimate objective in designing an efficient survey design is the 
maximization of accuracy per unit cost. To accomplish this, we need 
to know the cost structure of the survey. We can identify three types 
of costs in two-stage sampling: (i) overhead costs; (iz) costs that depend 
primarily on the number of Psus in the sample; and (ziz) costs that 
depend on the number of ssus in the sample. Since the overhead costs 
are fixed, they can be ignored in determining the sample sizes. The 
costs of sampling psus may consist of the costs of selecting, traveling 
to, locating each sampled psu, and setting up the measuring equip- 
ment. A simple cost function may be of the form 


mC, + mnCo, (4) 


where C, and C, are the costs of sampling a Psu and SSU, respectively, 
and mn is the total number of ssus sampled. Typically, however, the 
cost functions are more complex. In CosTUS, as we mentioned earlier, 
the cost of sampling a Psu varies from one PSU to another and depends 
primarily on the switching equipment in the central office. Further, 
the cost of sampling a telephone line (ssu) also depends on the 
switching equipment and so varies from one office to another. There 
is also a special cost structure in the transmission impairments survey 
in Example 2. Here, if trunks (edges) are selected by using two-stage 
sampling to determine the pair of end offices connected to the trunk, 
it is cost-efficient to select offices with many trunks rather than those 
with fewer trunks. 
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3.2 Determining the parameters in a simple situation 


Let us consider a simple situation to illustrate the concepts involved 
in determining the parameters of the optimum design. Suppose the 
number of ssus in each psu is the same and equals N, the cost function 
is given by eq. (4) and we use SRSWOR to select the units in both stages. 
We need to determine only m, the number of Psus to be sampled, and 
n, the number of ssus to be sampled from each selected psu. The 
variance of the H-T estimator can now be written (see Section 2.4) as 


2 i V3 
V2) =a-pZ+a-py~% (5) 


for some V, and V2. Here fi = m/M and fp = n/N. A comparison of eq. 
(5) with the cost function C = mC, + mnC; reveals that increases in m 
and 7 have opposite effects on the variance and costs. Also, it is clear 
that an increase in m results in greater reduction in the variance than 
a corresponding increase in 7. Since C, is typically much larger than 
C2, it is more costly to increase the size of the first stage sample than 
the size of the second stage sample. All of these factors must be taken 
into consideration in determining the optimum combination of m and 
n. 

As mentioned earlier, optimum levels of m and n can be determined 
by minimizing either (z) the variance subject to a cost constraint or 
(ii) the cost subject to some accuracy requirements. Both approaches 
yield essentially the same results. The problem can be formulated 
mathematically as minimizing a given function subject to a constraint. 
Standard numerical or analytical techniques (LaGrangian multipliers, 
Cauchy’s inequality) can be used to determine the optimum values of 
m and 7. In this particular simple situation, explicit expressions for m 
and 7 can be easily obtained. Suppose we want to minimize the cost 
subject to the condition that the variance eq. (5) does not exceed some 
value b. If we can ignore the finite population corrections in eq. (5), 
the optimum values of 7 and m can be obtained as 








Bop = Vo/Vi 
VC2/Ci 
and 
ry = LLL, 
where 


A = b/(Ci/Vi + Co/ V3)”. 
The total cost of the survey with these values of m and n is given by 
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Copt = b/A?. 


In practice, one should not be satisfied with just determining the 
optimum values of m and 7 without examining the behavior of the 
variance and cost functions near the optimum. Since preliminary 
estimates of costs and variances may only be approximate, the behav- 
ior of these functions in a neighborhood around the optimum should 
be examined. Relatively flat variance and cost functions near the 
optimum value indicate robustness against possible moderate errors in 
the input parameters. 


3.3 More general situations and the COSTUS example 


In most surveys, the situation is more complex than the one we have 
just discussed. For example, the psus will not necessarily be the same 
size and the cost function may be more complicated. Even in a general 
situation, the problem can be formulated in such a way that we can 
determine, either analytically or numerically, the optimum values of: 
m, the number of Psus to be selected; {II;}, the inclusion probabilities; 
and {n;}, the number of ssus to be sampled from each selected Psu. 
Some of these results for some special cost functions can be found in 
the literature.” 

We want to emphasize here the importance of simplifying the 
problem, whenever possible, by using reasonable approximations. In a 
complex situation where there are too many design parameters to be 
determined, it is difficult to appreciate the impact of unreliable input 
values. Reducing the number of parameters through the use of some 
practical guidelines usually provides us with a better understanding of 
the problem. We illustrate some of these ideas through the costus 
example. 

The PSUs in COSTUS are central offices and, as mentioned earlier, the 
cost of sampling the office and telephone lines (ssus) in the office 
depends on the type of switching equipment in the office. Since each 
office provides service in several classes, we have to sample lines from 
all the available classes in the selected offices. However, not every 
office provides service in every available class. Since we want to study 
the parameters of all the classes, we take the first-stage costs of 
sampling an office with service in only one class to be twice that of an 
office with service in two classes. Hence, the total costs of the survey 
can be written as 


m Ss 
TC=Y, {Du + > Duna, 
i=l C=1 
where D,; = D:/2;, Di; = the costs of sampling the ith office, 2; = 
number of classes in the ith office, D2; = costs of sampling a line from 
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the ith office, and nc; = number of lines to be selected from the ith 
office for class C (equals zero if office 1 does not have service in class 
C). Since the total cost of this survey is a random quantity, we 
minimize the expected total cost 


M _ s 
m b Zi (2. + Da ¥ na) | (6) 
i=] C=1 


where mZ; = II, the inclusion probabilities. 

We need to minimize eq. (6) subject to some accuracy constraints. 
In this study, the quantity to be estimated is the mean load (in ccs) 
during the busy hour, Yc. We require a relative error no larger than 
0.1 with probability 0.90 for each of the S classes, C = 1, ---, S. From 
Section 2.3, we note that this accuracy can be stated in terms of an 
upper bound on the mse of the estimator. We use the following 
approximate expression for the relative mse (rmse) to determine the 
design parameters: 





" _ [™M w%. /92, 
rmse(Yc) = Yo" | Wo (== 
i=1 Nci 


rar + (Yo - Por) | (7) 


The notation here is the same as in Section 2.2. The additional 
subscript indicates the class of service. This expression (which in fact 
equals the relative variance of the unbiased estimator) is the zeroth- 
order term in the Taylor series expansion for the rmse of the ratio 
estimator. By not taking into account the higher-order terms which 
include the correlation between the numerator and denominator of the 
ratio estimator, this expression, in general, overestimates the variabil- 
ity. However, it is simpler to use and the overestimation may be 
desirable in view of the unrealiability in preliminary estimates. 

Before determining the design parameters, we make two addi- 
tional simplifications: (i) replace SZ:/Y2 in the first component of eq. 
(7) by Vcz, a quantity that does not depend on the office 1; and (zt) 
replace (Yc: — Yc)?/Y@ in the second component by Va, also a 
quantity independent of the office. This is reasonable since a priori we 
do not expect much variation between these values and, in any event, 
we do not know each one of the individual values. (See also the 
discussion on prior information in Section 2.5.) 

Thus, we want to determine the design parameters which minimize 

M Wi, 
ies 


eq. (6) subject to 
(Vo + ve) < bec 
im mZ; Nci 


for some bc, C = 1, ---, S. Instead of determining m, {Z;} and {n;} 
from the optimality calculations, we only determine m and fic, the 
average number of ssus to be sampled from a selected Psu for each 
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class of service. Once m and 7ic are determined, we can allocate mnc, 
the total number of lines for class C, to each sampled office inversely 
in proportion to (Dz;)'””. We also select {Z;} in advance by taking them 
proportional to 


{N.:/Di;”"}, 
where 
Nui = Ye=-1 Nei. 


Once we substitute these values for {nc;} and {Z;} in the variance and 
cost functions, it is a relatively easy problem to find the values of m 
and nc that minimize the total expected cost subject to the accuracy 
constraints. Since there are only S + 1 design parameters involved, it 
is also easy to examine the behavior of the cost and variance functions 
near the optimum and investigate the sensitivity to errors in input 
values. 

When cosTus was implemented in a few jurisdictions, we also 
examined the advantage gained by using unequal probability selection 
schemes. Since we were using the conservative WR variance formulas 
for the unequal probability selection scheme, we found that the loss in 
“efficiency” from using SRSWoR of offices (with exact variance formu- 
las) was not substantial. This also simplified the computations consid- 
erably. 


IV. SAMPLING DESIGNS FOR POPULATIONS WITH PARTIAL VARIATE 
PATTERNS 


4.1 The problem of partial variate pattern (PVP) 


A multivariate population (for example, one with multiple charac- 
teristics) is said to exhibit a PvP if not all the variates can be observed 
from every unit in the population. In Costus, as we noted, not all the 
central offices provide service in every available class. In the survey of 
baseband transmission impairments in Example 2, not all carrier 
systems appear between each pair of central offices. It is easy to 
visualize many other studies, both within and outside the Bell System, 
where the populations exhibit pvp. The problem of PvP can be serious 
if there is great variation in the size of the universe corresponding to 
each variate. The usual sampling designs may not provide reasonable 
assurance that we can select a sample that will allow us to estimate 
the parameters corresponding to each variate with prescribed accu- 
racy. 

Let us consider some schemes for sampling in the presence of PvP 
(also see Ref. 13). Since the problem of PvP is present in one stage of 
the selection process only, we restrict our attention to sample selection 
in the first stage. Thus, suppose there are M units in the population, 
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of which Mc units have characteristic C, C= 1, --- , S. Let the sample 
size, determined by accuracy requirements, for characteristic C be mc. 
These sample sizes of course also depend on the particular sampling 
scheme used. 


4.2 Some sampling designs 
4.2.1 Modified simple multivariate sampling 


Let m = maxcmc and suppose we select a sample of m < M units, 
possibly using different selection probabilities for different units. This 
is the simple multivariate sampling scheme, intended for populations 
with no PvP. If mc denotes the number of sampled units with charac- 
teristic C, mc may be much smaller than mc and in some cases may 
even be zero. We can modify this scheme in a number of ways. Instead 
of selecting m = maxcmc units, we can select m* units, according to 
selection probabilities {II;}, where m* is determined so that the 
expected number of units in the sample is at least mc, C= 1, ---, S. 
This can be achieved by taking m* = maxcmc/pc, where pc is the total 
of the probabilities Z; = I1;/m for units with characteristic C. This can 
be justified if we view the selection of a unit with chracteristic C 
approximately as a binomial experiment with probability of success 
pc. This formulation can alternatively be used to determine m* such 
that, say 90 percent of the time, mc = mc, C= 1, ---,S. 


4.2.2 Combined multivariate sampling 


Here, we consider S universes, each universe corresponding to the 
units with characteristic C, C = 1, ---, S. We select an independent 
sample of size mc from each one of the S universes. We then observe 
every available characteristic from the units selected in all of the S 
samples. The total number of units selected in these S samples can 
vary between maxcmc and Y'2-1 mc. The main disadvantage of this 
scheme is that this number may be too large. However, we can exercise 
some control over this number. One possibility is to give higher 
selection probabilities to units with more characteristics than those 
with fewer characteristics (see Section 3.3). Alternatively, instead of 
selecting mc units from the universe corresponding to characteristic C, 
we can select a smaller number, mé, of units. This is because we expect 
to select some units, in addition to these mé units, with characteristic 
C from the remaining S — 1 samples. So, the number mé can be 
determined such that either on the average or with prescribed proba- 
bility, the total number of units with characteristic C exceeds mc, 
C = 1, ---, S. The binomial approximations discussed earlier can be 
used to determine the {mé}. 


1250 THE BELL SYSTEM TECHNICAL JOURNAL, SEPTEMBER 1981 


4.2.3 Stratified sampling 


We can also try to deal with pvp by stratifying the units so that, 
within each stratum, the units are internally homogeneous in some 
sense in terms of the pvp. We consider two stratification techniques 
here. 

In the first scheme, called variate stratification, the strata are 
determined in terms of the variates (characteristics). Suppose the 
variates are ordered so that the number of units with variate one is 
smallest, the number with variate two is next smallest, etc. Then, 
stratum one consists of all the units with variate one, stratum two 
consists of all units with variate two and not in stratum one, etc. If we 
now allocate the total sample size among the strata, we can estimate 
’ the parameters corresponding to all the variates, especially the “small” 
ones. However, this scheme is not foolproof in the sense that it is 
possible to construct examples where the selected sample does not 
contain any units with one of the variates. 

The second method, pattern stratification, is based on the variate 
pattern. Here, units with identical variate pattern, i.e., having the same 
set of characteristics, are grouped into a stratum. Unlike the variate 
stratification scheme, we can guarantee the required sample size for 
each variate in this scheme. However, this scheme suffers from the 
serious drawback that the total sample size may be too large, since the 
number of different strata (which is smaller than the sample size) can 
be as large as min(M, 2° — 1). 

In both these schemes, standard nonlinear programming techniques 
can be used to determine the sample size for each stratum to minimize 
cost subject to the variance constraints. 


4.2.4 Other methods 


It is possible to use sequential sampling schemes to ensure that we 
select a sample with a given number of units for each characteristic.* 
However, it is extremely difficult to determine analytically the selec- 
tion probabilities for most of these schemes. One simple sequential 
method that can be implemented is a two-stage simple multivariate 
sampling scheme in which a simple multivariate sample is supple- 
mented by a second-stage sample from the remaining units. Although 
the variance calculations become more involved, they are still tracta- 
ble. 

It also is plausible that ideas from the controlled selection method- 
ology can be applied to the selection of samples from populations with 
pvp.'*!>16 However, it is not clear how to characterize explicitly the set 
of all feasible samples here. Variance calculations also remain a difficult 
problem with controlled selection. 
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4.3 The design used in COSTUS 


The sampling design used in costus for handling Pvp will be 
described here. As the problem of PvP exists only in the first stage, we 
consider the selection of units in stage one only. 

While examining data from several jurisdictions for the different 
PVPs, we found that, in most cases, a class of service can be classified 
as either small or large in terms of the proportion of offices with service 
in that class. There were very few jurisdictions with medium-sized 
classes of service. 

Since the main concern in the presence of PvP is the ability to 
estimate parameters corresponding to the small classes of service, we 
decided to group all offices with services in these classes in stratum 1. 
A combined multivariate sampling scheme, which guarantees the 
required sample size from each class, is used to select offices from this 
stratum. Since the total number of offices sampled under this scheme 
may be large, we restrict the size of this stratum to be no larger than 
25 percent of the universe. 

We can use a simple multivariate sampling scheme to select a sample 
from the remaining offices. However, we first identify those classes 
with service in less than 50 percent of the remaining offices. The offices 
with service in these classes (and not in stratum 1) are grouped into 
stratum 2. The remaining offices are grouped into stratum 3. Simple 
multivariate sampling schemes are then used to select units in strata 
2 and 3. By doing this, we have reasonable assurance that the sample 
sizes for the classes that characterize stratum 2 are not too small 
compared to the required sizes. . 

Hence, we see that the sampling design for COSTUs is in fact a three- 
stage sampling design. In the first stage, the offices are grouped into 
three strata. Different sampling schemes are used in the different 
strata to select offices in the second stage. From each office selected in 
the second stage, telephone lines corresponding to each available class 
are selected in the third stage. 

The design we have used here for handling PvP incorporates specific 
features of some of the schemes discussed in Section 4.2. The stratifi- 
cation is based on considerations similar to those in the variate 
stratification scheme. It is, however, adaptive in the sense that it 
depends on the variate pattern in each universe. In our applications, 
we found that in many jurisdictions stratum 2 was empty and in some 
situations, where the problem of PVP is not serious, stratum 1 was 
empty. 

We arrived at the final design used in cosTuS by examining data 
from various jurisdictions for the different types of PvP to expect. This 
design, while not foolproof, provides a reasonable, practical solution to 
the problem at hand. 
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V. SAMPLING FROM NETWORKS 


In most surveys, we can treat the population under study merely as 
a collection of elementary units with no importance attributed to the 
interrelationships that exist among the units. In some situations, 
however, these relationships cannot be ignored and the selection of 
the sample is necessarily affected by the network of relationships that 
exist in the population. In this section, we briefly review some aspects 
of network sampling and discuss the sampling design used in Example 
2: 


5.1 Networks 


There are a wide range of surveys in the Bell System that deal with 
sampling from a network. Besides communication networks, network 
sampling also occurs in studies of other types of traffic flow and 
transportation facilities. A contact network or sociogram may repre- 
sent the interrelationships among a group of individuals, households, 
customers, etc. Other examples include similarity or dissimilarity struc- 
tures in cluster analysis and multidimensional scaling, where we want 
to compare a set of objects and group them into classes of similar 
objects. 

A network can be described in abstract terms with the aid of graph 
theory. An undirected graph (network) consists of a nonempty set V 
of elements called vertices (nodes) and a set of E of elements called 
edges. Each edge e of E is associated with a pair of vertices (i, 7). The 
edges may have several attributes associated with them. A network 
can also be represented by a matrix with the columns and rows 
representing the vertices. A one in the (i, 7)th cell of the matrix 
indicates that the vertices 7 and 7 are connected. In the survey of 
baseband transmission impairments discussed in Example 2, the ver- 
tices are central offices and the edges are trunks. In this case, there are 
many trunks and also different types of trunks between a pair of 
central offices. Several attributes, corresponding to the impairment 
characteristics, are associated with each trunk. 


5.2 Some sampling schemes 


The manner in which we have observational access to the elemen- 
tary units is the key to developing a reasonable sampling design. If we 
have a “frame” of all the edges in the graph from which we can select 
a sample of units, the problem is essentially one in traditional sampling 
theory. If no such frame is available and the structure of the relation- 
ship between the nodes must be discovered and explored during the 
course of data collection, the sampling design problem is quite differ- 
ent. Even in cases in which a complete listing of the edges is available, 
as in Example 2, cost considerations may dictate that a sample of 
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edges be selected by first sampling the nodes. Also, unlike traditional 
sampling where information about a unit can be obtained only by 
sampling and observing it, information about the relationship between 
several nodes may be obtained at any one of the nodes in network 
sampling. 

The field of sampling from networks has been considered by only a 
few authors so far.’’”” Most of the attention has been focused on 
surveys for which the structure of interrelationships is unknown and 
must be discovered. The references above deal mainly with estimating 
parameters that measure various aspects of these relationships. 

Goodman proposed the “snowball” sampling scheme for selecting 
edges (or pairs of connected nodes).”’ In this procedure, the survey 
proceeds from an initial sample of nodes by obtaining information 
about other nodes to which they are connected. The next step is to 
add to the sample some or all of these connected nodes, obtaining data 
from them as well as information about still other nodes to which they 
are connected. In an s-stage k-name snowball sample, this process is 
repeated for s stages and at each stage, k other nodes connected to a 
node already in the sample are selected. Goodman studies this scheme 
in detail under the assumption that the initial sample is selected 
through binomial sampling.” He also considers the case in which the 
k nodes are selected randomly at each stage. See also Ref. 1. 

To consider two other methods of network sampling, let us view the 
network as a matrix with the vertices corresponding to the columns 
and rows and the elements of the matrix corresponding to the edges. 
If we select a sample of nodes (rows/columns of the matrix), we can 
base our inference entirely on the sampled subnetwork that corre- 
sponds to the sampled rows and columns. This procedure (called 
subnetwork sampling) of selecting one or even several subsystems out 
of a number of subsystems is equivalent to traditional one-stage cluster 
sampling. It leaves open all questions about interrelationships between 
one cluster and another. In the partial network sampling scheme, we 
select a sample of nodes from the node set, and observe all the edges 
connected to one or more of the nodes in the sample. Estimation of 
the network characteristics using these two schemes is discussed in 
Refs. 1 and 18. 


5.3 Survey of baseband transmission impairments 


In this survey, there are a number of trunks of various types with 
each trunk associated with a pair of end offices. If we select a particular 
pair of end offices, it then becomes cheaper to select additional trunks 
from those trunks that terminate in either one of the two offices. This 
special cost structure implies that we need to select trunks (edges) by 
appropriately selecting offices (nodes) to which they are connected. 
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A multistage sampling scheme is used in this survey. A sample of 
primary offices, using probabilities proportional to some measure of 
size (the number of trunks), is selected in the first stage. A number of 
secondary offices are selected, again using probabilities proportional 
to some measure of size, from the set of offices connected to each of 
the primary offices. From every pair of end offices thus sampled, a 
number of trunks corresponding to each trunk type are selected using 
simple random sampling. The parameters of the sampling design (m, 
the number of primary offices, {m;}, the number of secondary offices 
and {ni}, the number of trunks of a particular type) can all be 
determined so that the total survey cost is minimized subject to some 
accuracy criterion. 

The two-stage sampling scheme used here to select the pair of end 
offices can also be viewed as a two-stage snowball sampling scheme. It 
is of course possible to use a k-stage snowball sample to select the 
offices. Optimality considerations relating to the number of stages and 
the sample size in a snowball sample have yet to be resolved. 


VI. SUMMARY 


We have reviewed various aspects of sampling from structured 
populations in this paper. The issues that have been selected for 
discussion, two-stage sampling from populations with multiple char- 
acteristics and sampling designs for populations with Pvp and network 
sampling, are common to many Bell System surveys. Thus, we hope 
that an exposition of some of the theoretical and practical considera- 
tions involved in dealing with these situations will serve other survey 
practitioners. Throughout the paper we have tried to balance theoret- 
ical considerations with practical guidelines gained from our own 
experience. Two recent Bell System surveys are used to illustrate the 
ideas disscused. 
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This paper describes a new scheme for allocating a data bus on a 
first-come-first-serve (FCFS) basis. When the devices connected to the 
bus fequest to become the bus-master, they are assigned distinct 
“ticket numbers” in the order in which the requests are generated, at 
which time they go into a wait state. When the bus is released by a 
device holding the ticket number n, it is then allocated to the device 
holding the ticket number n + 1. We discuss the conditions under 
which the scheme is a close approximation to the ideal FcFs scheme 
and evaluate its performance using simulation results. We also 
present two alternative hardware implementations of this scheme— 
one centralized and the other distributed. Because of its simple 
hardware implementation, the scheme is attractive for applications 
where a bus is shared, in an unbiased fashion, among a large number 
of devices. 


1. INTRODUCTION 


In computer systems, situations frequently arise where a resource is 
shared among several devices, but it can be used by only one device at 
a time. Scheduling such a resource to enforce mutual exclusion over 
its use is necessary if devices request the resource while it is being 
used or if the requests arrive simultaneously. A frequently encountered 
resource of this type is the data bus, which provides a communication 
path among the various devices connected to it. At any given time, 
there can be several devices receiving (or reading) information from 
the bus, but there can be only one device that has the privilege of 
transmitting information on it. Such a device is called the bus-master, 
and mutual exclusion among the devices wishing to become the bus- 
master is enforced by bus arbitration schemes. 

Devices requesting the bus while it is busy are made to wait until it 
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becomes available again. As soon as that happens, one among the 
waiting devices is allowed to become the bus-master. In most bus 
arbitration schemes, this choice is made without regard to the order in 
which the requests originally arrived; for example, daisy-chaining, 
device polling, and parallel priority resolution schemes. 

Some of the commonly used bus arbitration schemes have been 
reviewed by Chen and Thurber et al.” Among them, polling and daisy- 
chaining are most commonly used. Polling is suitable only for slow 
devices, because the waiting times from bus request to bus grant are 
quite long, as the devices can access the bus only during preassigned 
time intervals. Daisy-chaining is extensively used in several minicom- 
puters, such as the PDP-11s made by Digital Equipment Corporation 
(DEC).* Arbitration delay in this scheme may be quite long, since it is 
proportional to the number of devices connected to the bus. Further- 
more, by virtue of their location on the bus, the devices are assigned 
fixed priorities that are used for contention resolution. 

For faster bus arbitration, the recent computers designed by DEC 
and Honeywell, Inc. use distributed schemes.*” These schemes, and 
those described in Refs. 6, 7, and 8, use the same algorithm with 
different implementations. They are fast, modular, and flexible, but 
they, too, allocate fixed priorities to the devices connected to the bus. 

The major drawback of allocating fixed priorities to the devices is 
that the low priorities may have to wait indefinitely before being 
granted access to the bus if a few high-priority devices decide to use 
the bus frequently. They are effectivelly “locked out” from service. See 
Ref. 9 for a simulation-based quantitative analysis of these and other 
bus arbitration schemes. 

In this paper, we present a first-come-first-serve (FCFS) scheme that 
allocates the bus in an order that is a close approximation to that in 
which the devices request the bus. This scheme does not have the 
above-mentioned drawback of locking out a few devices from service, 
and it provides an equal grade of service to all devices. We first 
describe the scheme and then discuss two alternative hardware imple- 
mentations—one centralized and the other distributed. Before describ- 
ing the scheme, we briefly discuss the advantages of following the Fcrs 
allocation policy. 

Consider a data bus that is shared among several devices, and 
assume that (z) the devices request the bus with the same statistics, 
and (i) the bus is allocated for a fixed quantum of time for each 
request. The bus arbitration scheme should then have the following 
two properties: 


(i) It will minimize the idle time on the bus, so that the bus 
throughput is maximized. This is done by arbitrating for the 
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next bus-master concurrently with the bus usage and by reduc- 
ing the arbitration time, if it happens to be longer than the time 
quantum for which the bus is allocated. 

(zi) It will have the least disparity of service across requests and 
also across devices. The disparity of service across requests is 
represented by s, the standard deviation of the waiting times 
taken over all bus requests, and the disparity of service across 
devices is represented by S, the standard deviation of the 
average waiting times experienced by the individual devices. It 
is a simple matter to show that if an arbitration scheme does 
not prefer a device over any other, the average waiting times 
experienced by individual devices are all equal. Such schemes 
are called unbiased schemes, and S for them is zero. The ideal 
FCFs scheme is one such scheme. In the Appendix we show that 
the ideal rcrs scheme also attains the minimum value of s. 
Thus, under the assumptions stated above, the ideal FcFS 
scheme is a desirable scheme to be emulated in real systems. 


ll. DESCRIPTION OF THE SCHEME 


Let there be N devices connected to the bus. In order to ensure that 
the devices gain bus control one at a time and in the order in which 
they requested it, we propose a scheme that is very similar in essence 
to that used in many supermarkets. As customers walk in, they pull 
out a numbered ticket from a machine that dispenses sequentially 
numbered tickets. When the server becomes free, he or she waits on 
the customer with the ticket one number higher than that of the last 
customer served, thus, providing equitable service to all customers. 

Our scheme is based upon two essential pieces of information: “next 
number to be served” (NNS) and “next number available” (NNA). This 
information can be maintained in a centralized or distributed fashion, 
as we discuss in Section III. In addition, each device has a register, 
called the ticket register, to store the ticket number assigned to it 
when it requests bus mastership. How these ticket numbers are as- 
signed is discussed later; let us first see how they are used. As soon as 
the bus is available, each device compares its ticket number with the 
NNS, and the device that finds the match becomes the bus-master. As 
we explain in the following discussion, there can be only one device 
whose ticket number matches NNS. Sometime before the bus is avail- 
able again, NNS is incremented by one. This incrementing is done 
modulo NTICKETS, so that the ticket numbers range from 0 to 
(NTICKETS-1). To ensure that devices have distinct ticket numbers, 
we must have NTICKETS = N, where N denotes the number of 
devices. 
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Now we consider how the ticket numbers are assigned. This is done 
using NNA. When a device wants to become the bus-master, it copies 
the NNA into its ticket register. Then, NNA is incremented by one 
modulo NTICKETS, thus, ensuring that sequentially increasing ticket 
numbers are “dispensed” in the range from 0 to (NTICKETS-1). Of 
course, while the copying and incrementing operations are being done, 
no other device should be allowed to copy the NNA. If the NNA is 
copied, either there will be two devices with the same ticket number, 
and confusion will ensue as they both will become bus-masters at a 
later time, or there will be a device with an invalid ticket number that 
is outside the above-specified range, and that device would never be 
able to access the bus, as NNS will never be equal to the invalid ticket 
number. The accesses to NNA to receive ticket numbers should, there- 
fore, be mutually exclusive. 

Thus, in our ticket assignment scheme, achieving mutual exclusion 
for the bus depends on achieving mutual exclusion at a lower level— 
that of NNA. The second mutual exclusion is achieved by using one of 
the existing arbitration schemes; for example, simple daisy-chaining 
(spc), rotating daisy-chaining (RDC), modified device polling (MDP), 
dynamic parallel priority resolution (DPPR), etc. See Ref. 9 for a 
detailed description and comparison of various bus arbitration 
schemes. 

The duration for which NNA is allocated to a device is the time it 
takes to copy NNA into its ticket register. This duration is very short— 
typically, a few gate delays. Thus, the time required to assign a ticket 
number is essentially the time spent in arbitrating for the use of NNA. 
Whenever this time is short, as compared to the time for which the 
main bus is allocated, our scheme would be a close approximation to 
the ideal rcrs scheme. This is because the devices that request the 
bus while it is busy are quickly assigned ticket numbers and put into 
a waiting state. Thus, the scheme remembers the order in which the 
requests arrive. 

It is also possible that some devices will request the bus while NNA 
arbitration is in progress or while NNA is in use. In such cases, 
depending upon the NNA arbitration scheme, one among these devices 
is allowed to copy the next NNA, and they may or may not receive the 
ticket numbers in the temporal arrival order of their requests. Thus, 
for the overall scheme to be unbiased, the NNA arbitration scheme 
must treat the devices in an unbiased way. This is desirable because it 
is a necessary condition for making the overall scheme a close approx- 
imation to the ideal Fcrs scheme. 

To summarize, the NNA arbitration scheme should be fast and 
unbiased. We use the criteria in choosing the NNA arbitration scheme. 
In addition, to judge how close the overall scheme is to the ideal Fcrs 
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scheme, we use s and S. The smaller the value of s, the better the 
approximation, since s is minimum in the ideal case. Similarly, the 
smaller the value of S, the better the approximation, since S is zero in 
the ideal case. 

We now examine the SDC, RDC, MDP, and DPPR schemes mentioned 
earlier with regard to their desirability as NNA arbitration schemes. 

In spc, the central arbiter sends out a daisy-chained NNA-grant 
signal. If a device does not want to access the NNA, it lets the signal 
pass through; otherwise, it stops the signal and then accesses NNA. 
Thus, the devices closer to the arbiter are preferred over those farther 
away from it. This tends to make the overall scheme, named FCFS/ 
SDC, a poorer approximation to the ideal Fcrs scheme. 

In RDC, the device that accessed NNA last acts as the arbiter for the 
next arbitration cycle and sends out the NNA-grant signal. On the 
average, all the devices are given equal treatment, and the average 
arbitration time is the same as that for spc. Therefore, the overall 
scheme, FCFS/RDC, is expected to be a better approximation of the 
ideal FcFs scheme than the FCFS/SDC. 

In MDP, there is no central arbiter, and the daisy-chained NNA-grant 
signal keeps travelling from device to device in a cyclical fashion. 
Devices wishing to access the NNA wait for the grant signal to arrive, 
stop the grant signal temporarily, access the NNA, and then release the 
grant signal. All the devices receive unbiased treatment. The perform- 
ance of FCFS/MDpP is, therefore, expected to be similar to that of FcFs/ 
RDC, and their hardware implementation is also quite similar. 

In DPPR, devices are assigned priorities which change after each NNA 
arbitration cycle. As the arbitration starts, all the devices that need to 
access the NNA put their priorities on a common priority bus. Then, 
each device removes itself from the contention if its priority is lower 
than the composite priority on the priority bus. This eliminates all but 
the highest priority device, which then accesses the NNA.” The dynamic 
assignment of priorities in this scheme can be done in a variety of 
ways, but here we assume it is done so that the order in which devices 
win arbitration is essentially the same as that of RDC (the same priority 
assignments emulate MDP also). Initially, the ith device is given the 
priority i, and after each arbitration, priorities are cyclically rotated so 
that the device that won the last arbitration gets the priority one. All 
the devices are treated equally; however, the average arbitration time 
for DPPR is much smaller than that of RDC or MDP, because there is no 
daisy-chained signal involved that gets delayed while passing through 
each device (by as much as 4 gate delays per device). Thus, FCFS/DPPR 
is a better approximation to the ideal FcFs scheme than FCFS/RDC and 
FCFS/MDP. The disadvantage is that FCFS/DPPR requires more hard- 
ware than FCFS/RDC or FCFS/MDP. 
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The conclusions drawn above are supported by the values of s and 
S obtained through simulation, which are shown in Tables I and II, 
respectively. (These statistics have been borrowed from Bain and 
Ahuja.’) Both tables include two cases: (i) when the schemes discussed 
above are used for NNA arbitration in the ticket assignment scheme, 
and (ii) when they are used to arbitrate for the main bus itself. Note 
that the values of s and S for ticket assignment schemes (column 1) 
are smaller than for the others (column 2); therefore, they are better 
approximations to the ideal rcrs scheme. Similarly, among the ticket 
assignment schemes, FCFS/DPPR is the closest approximation to the 
ideal FcFs scheme. Through simulations, we also observed that as the 
number of devices is increased, the performance of FCFS/DPPR rapidly 
converges to that of the ideal Fcrs scheme, but the performance of 
other schemes diverges significantly from that of the ideal FcFs 
scheme. Therefore, FCFS/DPPR is an attractive scheme when a large 
number of devices (approximately 16 or more) share a common bus. 


Il. IMPLEMENTATION OF THE TICKET ASSIGNMENT SCHEME 


In this section, we describe and compare two implementations of 
the ticket assignment scheme. In the first, the NNA and NNS are 
centralized, and in the second, they are distributed. We consider only 
the FCFS/MDP scheme, since the implementations with different NNA 
arbitration schemes are quite similar. 

Figure 1a shows an implementation of FCFS/MDP in which the NNA 
and NNS counters are centralized, and the devices access them through 
the NNA and NNS buses. If a device requests access to the main bus, it 
waits for NNA-GT, the cyclically daisy-chained NNA grant signal, to 


Table |—Table of s, the standard 
deviation of the weighting times taken 
over all requests to the bus. X denotes the 
schemes in the first column and FcFs/X 
denotes the ticket assignment using X for 
NNA arbitration. 


s for FcFS/X s for X 

Xx (us) (us) 

SDC 19.32 30.26 
RDC 1.476 3.214 
MDP 1.512 3.230 
DPPR 1.368 3.159 
Ideal FcFS 1.112 


Note: Simulations were carried out for 32 inde- 
pendent devices, each device requesting the bus with 
uniformly distributed interrequest times between 0.4 
and 19.6 ps, with the average interrequest time of 10 
ps. The bus was allocated for 0.4 us for each request. 
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Table IIl—Table of S, the standard 

deviation of the average weighting 

times experienced by the individual 
devices. Simulations were carried out 
under the same conditions as shown 


in Table I. 
S for FcFs/X S for X 
Xx (ns) (ns) 
SDC 524000* 132600t 
RDC 68.0 199.9 
MDP 63.7 155.8 
DPPR 45.0 153.9 


* Devices 26 through 32 did not get service. 
+ Devices 23 through 32 did not get service. 


arrive. The device then holds the NNA grant signal, transfers the data 
on the NNA bus to its ticket register, signals on the INC-NNA line to 
increment the NNA counter, and then releases the NNA grant signal. 
(See Fig. 1b for a detailed circuit diagram.) After the device has a 
ticket number, it waits until the contents of its ticket register are the 
same as the data on the NNs bus and the bus busy line, BB, is negated. 
When that occurs, it asserts BB, becomes the bus-master, and signals 
on the INC-NNS line to increment the NNs counter. After using the bus, 
it simply negates the BB line. The BB line permits incrementation of 
NNS to proceed while the main bus is being used. 

Figure 2 shows a distributed implementation of the above scheme in 
which each device has its own NNA and NNS counters. The NNA and 
NNS buses are eliminated, and the INC-NNA and INC-NNS lines are used 
to keep the various NNA and NNS counters up to date. The initial 
values of all the ticket registers, NNA counters, and NNS counters are 
0, 1, and 1, respectively. 

When a device wants a ticket number, it executes the following 
sequence of steps: 


(i) Waits for NNA-GT to arrive, and captures it on arrival. 

(ii) Initiates steps (iii), (iv), and (v) when INC-NNA becomes false, 
and does nothing before then. 

(tit) Shifts the contents of NNA counter into the ticket register. 

(tv) Signals all the other devices to increment their NNA counters 
by asserting INC-NNA line. The device also signals itself to do 
the same. The INC-NNA line is negated after a long enough time 
to allow the NNA counter to finish incrementing. 

(v) Releases the NNA-GT signal. 


Notice that only steps (i) and (v) depend on the NNA arbitration 
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Fig. la—An implementation schematic for FrcFs/MDP where NNA and NNS are cen- 
tralized. INC-NNA and INC-NNS are used by the devices to increment NNA and NNS, 
respectively. NNA-GT is the cyclically daisy-chained grant signal for accessing NNA. 


scheme being MDP; they are replaced by a different set of steps for 
different NNA arbitration schemes. All other steps, including those 
given below, are independent of the NNA arbitration scheme used. 
When a device receives the INC-NNA signal, it simply increments the 
NNA counter. 

In the distributed implementation, gaining control of the main bus 
is similar to that in the centralized implementation: 


(t) After receiving the ticket number, wait until the contents of 
the NNS counter and the ticket register are the same, and INC- 
NNS and BB are false. When that occurs, initiate steps (iz) 
through (uv); do nothing before then. 

(tt) Set BB to true. 

(iii) Signal all other devices to increment their NNS counters by 
asserting the INC-NNS line. The device also signals itself to do 
the same. The INC-NNS line is negated after a long enough time 
to allow the NNSs counter to finish incrementing. 

(tv) Use the main bus. 

(v) Release the bus by setting BB to false. 


As a device receives INC-NNS, it increments its NNS counter. The 
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detailed circuit diagram for this is similar to Figure 1b, except that 
each device has NNS and NNA counters of its own. 

The distributed implementation has two advantages over the cen- 
tralized implementation. 
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Fig. lb—The detailed circuit diagram associated with the schematic Fig. la. The 
circuit enclosed within the broken lines is contained in each device. F/F denotes flip- 
flop. The bus-request F/F and NNA-GT F/F capture the grant signal. The one-shot A 
generates the outgoing grant signal, the negative of which is also used to signal on the 
INC-NNA line. The one-shot B generates the INC-NNA signal. Notice that the buses use 
negative logic. 
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Fig. 2—An implementation schematic for FCFS/MDP where NNA and NNS are distrib- 
uted. Each device has NNA and NNS registers. INC-NNA and INC-NNS signals are used to 
keep all the NNA and NNS registers up to date. The NNA and NNS buses have been 
eliminated. 


(t) It has fewer lines, since it does not need the NNA and NNS buses. 

(ti) Devices do not have to wait for voltage levels on the NNA and 

NNS buses to settle down, as NNA and NNS are available from 

their local counters. The longer the bus, the more significant this 

advantage because the settling time of voltages on the bus is propor- 
tional to the length of the bus. 


The distributed implementation has two disadvantages as compared 
to the centralized implementation: 


(t) In order to introduce new devices in the system, their NNA and 
NNS counters must be current with those in other devices. 
Although not always satisfactory, this can be done by stopping 
the system momentarily to reset the counters. 

(tt) The scheme will malfunction if any one of the counters mal- 
functions. Depending upon the reliability of the hardware, this 
disadvantage may not be serious. 


Thus, since neither implementation is unequivocally superior to the 
other, the final choice should be made depending upon the require- 
ments of the application at hand. 


IV. SUMMARY 


We presented a Fcrs bus arbitration scheme that is based upon 
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assigning ticket numbers to the devices as they request the bus. The 
arbitration for the main bus essentially depends upon the arbitration 
for the next available ticket number. Several schemes for the latter 
arbitration were considered, and their impact on the overall scheme 
was examined using the standard deviation of wait times of all requests 
and the standard deviation of the average weight times of devices. 
Using simulation results, we showed that the overall scheme is the 
closest approximation to the ideal Fcrs scheme, when the lower level 
arbitration is performed by the dynamic, parallel priority-resolution 
scheme; the resulting overall scheme is called FCFS/DPPR. Two alter- 
native implementations, one centralized and the other distributed, of 
the overall scheme were described. 
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APPENDIX 


In the following, we show (i) that the ideal rcrs has the minimum 
value of s, the standard duration of waiting times of all the requests, 
and (iz) that all disciplines of serving the requests have the same, 1, 
the average waiting time of requests. 

Consider the idealized arrangement where the incoming requests are 
put in a queue in their order of arrival, and the server always picks the 
first—the oldest—element, in the queue. This is the ideal Fcrs scheme. 
If at any time, the elements in the queue are permuted, we obtain 
deviations from the ideal case. 

Let w; be the waiting time of the zth request when the queue is not 
disturbed. Then, for the ideal Fcrs scheme, 


=. of 
w= iy 
nue 
and 
i 
2 —\2 
s°=—)' (wi — w)’,; 
= (wi @) 
where N is the total number of requests. 
Since any permutation can be expressed as a composition of a 
number of permutations that exchange two elements, we show that 
the value of w remains the same and that the value of s” is increased, 


if two elements in the queue are interchanged. For simplicity, we 
assume that the 7th and the (i — 1)st elements are interchanged. A 
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similar argument holds for the general case also. The new waiting 
times for these two elements are 


wi = uw; — t, 
and 
Wi-1 = Wi-1 +t, 


where ¢ is the service time for each request. Hence, the difference 
between the new and the old values of w is 


3. 1 
Aw = NW (wi-1 + wi) - Ni (wi-1 + wi) 


= 0. 


Also, the difference between the new and the old values of s? is 
1 
As? = W [(wi-1 — w — Aw)? + (wi — © — Aw)? ] 


= = [(wi-1 — W)? + (w; — w)*] 


Beane — wW;) 
N w—1 tse 


Note that the maximum value of w; occurs when the 7th request arrives 
in the queue immediately after the (z — 1)th request. If the ith request 
comes later, then the server services some requests in the meantime, 
thus, reducing the waiting time of the th request. Hence, 


W; S wWi-1 + . 
This gives us 
As’? =0, 


where the equality occurs only when the ith and the (i — 1)th requests 
come at the same time. Hence, the ideal rcrs scheme has the minimum 
value of s”. 
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Totempole and MOS Gates and Tristate 
Devices 
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The two logic values, 0, 1, and the unknown, are not sufficient for 
accurately simulating the behavior of TTL totempole and MoS gates 
and tristate devices. Furthermore, the classical fault modes (output 
stuck and input open) are not sufficient to cover the faulty behavior 
of Mos devices. A previous solution to the simulation modeling re- 
quired the addition of pseudo gates, which have no physical meaning. 
This paper develops methods of modeling fault-free and faulty tristate 
devices for logic simulation. The model does not require any addi- 
tional circuitry, but the existence of a simulator capable of simulating 
any number of logic values is assumed. 


1. INTRODUCTION 


A component finding wide usage in the bus-oriented architecture of 
today’s computer systems is the tristate driver. A typical arrangement 
is shown in Fig. 1. In this arrangement there are several drivers. Only 
one driver can be enabled at a time and “talk” to the bus. The receivers 
capture the information on the bus. 

In transistor-transistor logic (TTL) technology, tristate devices allow 
bus wiring, previously obtained only with conventional open collector 
output TTL. They also allow the use of active pull up to charge the 
large capacitances associated with the bus. This feature, not available 
with conventional bus wiring technique, speeds up the operation of 
the bus. In Mos technology, similar effects can be achieved with lower 
power requirements. Here, we shall consider TTL totempole, cmos, 
PMOS, and NMOS bistate and tristate devices and show the similarities 
and differences. 

In tristate technology, several malfunctions due to the presence of 
faults or to a wrong utilization of the bus may occur. These malfunc- 
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Fig. 1—A tristate bus system. 


tions may invalidate test results or damage the components. Therefore, 
it is important to simulate accurately the operation of faulty and fault- 
free circuits containing buses, and other tristate devices. 

It has been shown that the classical fault modes (output stuck, input 
open) are not sufficient to cover the faulty and fault-free behavior of 
cmos devices.’” One attempt has been made to map these fault effects 
into classical stuck-type faults by adding circuitry to the fault-free 
circuit. This additional circuitry is used to provide the faulty and fault- 
free circuits with memory properties, which exist in CmMos devices 
under certain conditions. This mapping allows the use of a fault 
simulator, which simulates only classical faults, for simulating faults in 
CMOS devices. In fact, the limitations of the available simulator was a 
constraint on the proposed modeling. Although modeling of the fault- 
free tristate devices also used similar added circuitry, the effect known 
as overlap or bus contention which may damage the devices, was not 
covered by this model. 

This paper develops methods of modeling fault-free and faulty 
tristate devices for logic simulation, without additional circuitry. How- 
ever, it assumes the existence of a simulator capable of simulating any 
number of logic values. Both the memory properties of Mos devices 
and the effects of bus contention are shown to be modeled accurately 
by the proposed method. 


1.1 TTL tristate technology 

Consider the tristate inverter of Fig. 2, in which A is the data input 
and E£ is the enable lead. The device is enabled and acts as an inverter, 
when £ = 0 as shown in Table I. When the inverter is disabled, it 
assumes a high impedance, namely the impedance between V, and 
Vcc, and the impedance between V, and Vg is extremely large. V., Vcc, 
and V¢ are the output, supply and ground voltages, respectively. 
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Table |—Truth table for tristate 





inverter 
A E Vo 
0 0 1 
0 1 Disabled 
1 0 0 
] 1 Disabled 
E 
A Vo 
Fig. 2—Tristate inverter. Fig. 3—TT tristate inverter model. 


The TTL tristate inverter can be modeled as the connection of two 
functions T; and T; (Fig. 3) which are controlled by lines A and E. T; 
and T. can be either conducting (on) or nonconducting (off) and they 
operate according to Table II. The device is in the high-impedance 
state when both T; and T» are off. In a tristate bus system, several 
tristate devices are wired together and the system operates safely if at 
most one device is enabled at one time (Fig. 1). 

Two problems have emerged in tristate bus technology and they are 
associated with the structure of the tristate devices. The first difficulty 
concerns a disabled tristate device and its ability to source or sink 
current depending on the value of its output voltage. These two cases 
are illustrated in Fig. 4. One can identify a voltage Vis, such that, if V. 
> Vin, then T2 acts as a current source, and if V, < Viz, Ti acts as a 
sink. Threshold voltage V1, is a voltage between Vcc and Vg, which is 
determined by the output properties of the device. 


Table Il—Truth table for inverter model 


A E T; T2 
0 0 Off On 
0 1 Off Off 
1 0 On Off 
1 1 Off Off 


If a receiver is present on a bus—i.e., it is connected to V,—the 
driver will source or sink current and, as a result, the output V, may 
reach the input threshold voltage V,, of the receiver. The output of a 
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(a) (b) 
Fig. 4—Disabled bus. (a) Sink. (b) Source. 


receiver with the input voltage equal to Vi, is unknown and in fact the 
output may oscillate due to small variations around Viz. Normally, 
after all the driving devices become disabled, the existence of the 
leakage current J; will destroy the previous logic value of a bus, and 
the bus will “float.” 

A second problem occurs when at least two tristate devices feeding 
a bus are simultaneously enabled and are in opposite active logic states 
(Fig. 5). Under this condition, the bus voltage may be anywhere in the 
range between the active logic levels and the currents may become 
extremely large (Fig. 5b). The actual value of V, and J, can be 


INVERTER Dy, INVERTER Dg 





(a) (b) 
Fig. 5—Bus conflict. 
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determined from the current-voltage characteristics of the two devices. 
This condition, called overlap, may cause excessive device heating 
resulting in device failure or slowly degrade the device, causing a 
decrease in life. 

In simulation, it is important to correctly model the effects of these 
two problems, so that a simulation user can be warned of the existence 
of potential difficulties. For instance, a primary input-output bus 
should be disabled (floating) before a test can be applied to it and 
enabled before the result of a test can be read from it. Also, an incorrect 
design or test sequence may cause overlaps on buses and the simulation 
should produce a warning. 


ll. TTL TOTEMPOLE, CMOS GATES, AND TRISTATE TRANSMISSION 
DEVICES 


2.1 Pull-up and pull-down functions 


Consider the cmos and TTL implementations of an inverter (Fig. 6). 
They have a common structure, which can be generalized by the 
diagram of Fig. 7 for multiple input gates. This structure is composed 
of a pull-up function (PUF), a pull-down function (PDF), and an inte- 
grator (I). The pUF and PDF depend upon the input values x1, +++ , Xn 
and the integrator produces the output Y depending upon Yy and Yp. 
Also the PUF and the ppF are complementary and cannot produce the 





Voc 
Vop 
PULL UP PULL UP 
P 
x 
Vo 
x Vo 
N 
PULL DOWN PULL DOWN 
Vss Vg 


(a) (b) 


Fig. 6—(a) CMOs inverter. (b) TTL totem pole inverter. 


LOGIC SIMULATION MODELS 1275 





Fig. 7—General model for cmos and TTL totem pole devices. 


same logic value under normal circumstances. By convention, Yu or 
Yp has the value 1(0), when the PUF or the PDF is on(off). The 
integrator IJ has the behavior of Table III, except in the case of 
malfunctions. As an example, consider the CMOS, NAND gate of Fig. 8. 


The junctions P; and P>2 realize a NAND function, whereas N; and Ne 
realize an AND function. 


Table Ill—Truth table for 


integrator 
Yu Yop Y 
0 0 Impossible 
0 1 0 
1 0 1 
1 1 Impossible 


2.2 Tristate devices 


A general CMOS or TTL tristate device can be modeled as in Fig. 9. 
The symbol E represents an enable line, and both PUF and PDF can be 
simultaneously disabled. Under normal conditions, the PUF and PDF 


cannot be simultaneously active. The integrator I is described in Table 
IV. . 


Table |!V—tTruth table for tristate 


integrator 
Yu Yb Y 
0 0 High impedance 
0 1 Logic 0 
1 0 Logic 1 
1 1 Impossible 
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Usually, tristate devices are used in the mode shown in Fig. 10. In 
the illustration, FE, Eo, --- E,, are the enabling lines. There are several 
interesting cases, namely 

(i) All the devices are disabled. 

(tt) One device is enabled. 
(iii) Two or more devices are enabled with opposite logic values. 
These three cases lead to different impedance situations on the bus 
and they are represented in Table V. In the context of simulation, it is 
possible to find a fourth situation, when the impedances are not known. 
This can be caused by an unknown value on the enable line E of a 


device. 
Vpp = "1" 
x 
x4 P, P2 
N2 
Ny 


Vg = “0” 


Fig. 8—CMOS NAND gate. 





Fig. 9—General model for tristate devices. 
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Table V—Device states and bus impedances 


Devices Impedance 
All disabled High pur, high PpF 
One enabled High Pur (PDE), low PDF (PUF) 
Two enabled (conflict) Low PUF, low PDF 


DRIVERS 





Fig. 10—Tristate devices connected to bus. 


2.3 CMOS dynamic properties 


The main difference between CMOS and TTL tristate devices is that 
the leakage currents in CMOS are extremely small compared to TTL 
leakage currents. Input currents for cMos are also small. Therefore, if 
a CMOS device is first enabled and then disabled, the small capacitances 
on the buses will remain charged for a long period of time and the bus 
will appear to receivers as if it were remaining at the same logic value. 
Ultimately, the capacitance will be discharged, but if the rate of 
operation is sufficiently fast, the discharge time can be considered as 
infinite and the bus displays memory properties. However, if forced to 
an active logic value, the bus will immediately reach this value inde- 
pendent of this charged capacitance. In TTL tristate devices, the 
leakage currents being large, the discharge time becomes small and no 
memory is displayed. 


2.4 Logic values and impedance 


It should be clear from our preceding discussion that accurate 
simulation modeling of tristate devices requires two distinct concepts: 
logic value and impedance. Three logic values, 0, 1, and u, are widely 
used in simulation, the symbol u being used to represent unknown 
signal values.? Unknown signal values may be present because the 
initial values of some leads may be unknown or because of races or 
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oscillations. The effects of impedance on circuit behavior depends on 
the technology. For instance, a high impedance may appear as an 
unknown logic value in TTL tristate technology. On the other hand, an 
output which has a high impedance in cmos technology will remember 
the logic value before the gate was disabled. Similarly, a conflict on a 
bus may appear as a 0, 1, or u depending on the technology. 

During simulation of circuits containing tristate devices, it is impor- 
tant to be able to detect special situations like bus conflict. Tests that 
cause bus conflicts may result in damage to the devices and must be 
avoided. In the tester environment, the state of an output bus in the 
high-impedance state may be altered by the tester, invalidating the 
test. These considerations lead to the representation of the state of a 
line by a pair composed of the impedance value and the logic value. 
There will be four possible impedance values (Table VI). Therefore, 
we obtain 12 combinations of impedance and logic value (‘Table VII). 


Table Vi—Impedance values 


Impedance 
Representa- 
PUF/PDF Impedance tion 
Both off High H 
One on, one off Regular R 
Both on Conflict C 
One or both unknown Unknown U 


Table Vil—Combinations of impedance and logic values 


Pair Description 
1 R/0 Logic 0 
2 R/1 Logic 1 
3 R/u Unknown with a low impedance 
4 H/0 High impedance with previous state memory 
5 H/1 High impedance with previous state memory 
6 H/u High impedance with unknown previous state memory 
7 C/0 Conflict with logic 0 effect 
8 C/1 Conflict with logic 1 effect 
9 C/u Conflict with logic u effect 
10 U/0 Unknown impedance, logic 0 effect 
ll U/1 Unknown impedance, logic 1 effect 
12 U/u Unknown impedance, logic u effect 


These 12 logic combinations, which we shall call logic values in the 
context of simulation, represent a detailed analysis of tristate devices 
and any one of these corresponds to a possible situation. Two cases 
are illustrated in Fig. 11 and both cases display memory properties. In 
the first case (Fig. 1la), the enable line goes from 0 to u and the output 
goes from R/O to U/0 (unknown impedance). In the second case (Fig. 
11b), the enable line goes to 1 and the impedance goes from R to H, 
with the same logic value. 
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E=0-u E=0—>1 










Y= R/O U/0 Y=R/1—-H/1 


(a) (b) 


Fig. 11—Determination of impedance. 


It is possible to reduce the number of values in Table VII at the 
expense of some information. The result is Table VIII, which shows 
two possible sets of logic values for TTL tristate devices. Z and a are 
synonyms for the pairs H/u and C/u, respectively. In set 2, pairs 3 and 
12 are differentiated. Pair 12 is called a potential conflict (a*) and can 
occur in various situations (Fig. 12). In this case, the simulation will 
declare a potential bus overlap. In set 1, pairs 3 and 12 are not 
differentiated and some information may be lost. In TTL technology, 
the unused combinations correspond to impossible situations (e.g., pair 
10) or to unpredictable situations (e.g., pair 8). 


BUS = a* 





(a) (b) 


Fig. 12—Potential bus conflict. 


Table Vill—Two sets of logic values for 
TTL tristate devices 


Number Pair Set 1 Set 2 
1 R/0 0 0 
2 R/1 1 1 
3 R/u u u 
4 H/0 Unused Unused 
5 H/1 Unused Unused 
6 H/u Z Z 
7 C/0 Unused Unused 
8 C/1 Unused Unused 
9 C/u a a 
10 U/0 Unused Unused 
11 U/1 Unused Unused 
12 U/u u a* 
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In the case of cmos technology, several additional values become 
meaningful (‘Table IX). The value ZO (Z1) is used when a driver having 
the value 0(1) is disabled and remembers the previous logic value. The 
value u0(ul1) is used when a driver was producing a 0(1) on a bus and 
the enable line becomes unknown. In all the three sets, pairs 7 and 8 
could be used if the actual conflict voltage can be positioned as a 0 or 
a 1 when added structural knowledge is available. 


Table IX—Logic values 
for CMOS devices 
Set 3 


ae) 
E. 


a 
NOR OW ONDA WONH 

C 

=] 

i=] 

n 

v0) 

Q 


The sets of values given in Tables VIII and IX can be used instead 
of impedance/logic value pairs. They are more economical in computer 
storage and more general; on the other hand, the pair representation 
may be more efficient, since the impedance is ignored in most of the 
gate evaluations (except the bus). 

We shall illustrate the use of the pairs for a CMOs driver-inverter 
(Fig. 2). The behavior of the inverter is represented in Table X. A and 
E are the input and enable lines, respectively, and Y is the output of 
the device. The symbol x represents a “don’t care” value. 


Table X—Impedance-logic value table for CMos- 
driver inverter 


Previous 
Valueof Y AE=x0 AE =01 AE = 11 AE = uu 
0 H/0 R/1 R/0 U/u 
1 H/1 R/1 R/0 U/u 
u H/u R/1 R/0 U/u 


A bus with any number of drivers may be calculated iteratively 
using Table XI. This table is symmetric with respect to the main 
diagonal. Its construction is illustrated by the case of three inverters 
producing R/0, R/1, and H/1, respectively, and wired to a bus. The 
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first pair produces C/u, and C/u is combined with H/1 to produce 
H/*, which can be approximated by H/u. 


Table Xl—Impedance-logic value table for tristate bus 





* Unknown (u) or technology-dependent value. 


2.5 Refinement of unknown impedance values 


Given that there are three basic impedance values, R, H, and C, the 
possibility of the enable signal being unknown during simulation 
introduces indeterminacy in the simulated impedance values. Seven 
impedance values can be used to represent the three known values 
and the four cases, where the impedance cannot be uniquely deter- 
mined. The seven values are: 


h=H 

Ih=R 

Ir=C 

i= HorkR 

Is = H or C 
Ig=RorC 
I;=HorRorC 


Figure 13 shows how the indeterminate simulated impedance may be 
generated. The impedance values are obtained by computing the 
impedances for a combination of the unknown signal values. For 
example, in Fig. 13c, it was possible to obtain J; because it was known 
that £; = E,=0or 1. 

These seven impedance values preserve some information that 
would otherwise be lost. For instance Js, Ig, and J; represent a potential 
overlap, whereas J, is definitely not an overlap. However, the overhead 
of dealing with a multiplicity of pairs may not be justified by the gain 
of information (21 impedance/logic value combinations). 
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uN 
Cc 


Ey =u E, Fu E 
R/1 or H/u 


! 


Eg =u Ig/u Eo=u Iyz/u Tg /u 


R/O 


(a) (b) (c) 


Fig. 13—Generation of indeterminate impedances. 


In the above analysis, indeterminacy in the impedance and logic 
value are treated separately. This also results in some loss of infor- 
mation, which becomes apparent from the computed output of the 
upper tristate inverter in Fig. 13a. Although its output is known to be 
R/1 or H/u, it will be represented by J,/u. Eliminating this problem 
would require creating one logic value for each subset of the set: {R/ 
0,R/1,R/u,C/0,C/1,C/u,H/0,H/1,H/u}, excluding the empty subset. 
This system would have 2°—1 or 511 logic values, which is impractical. 
In Fig. 14, the results would then become {C/u,R/0} for case a, {C/ 
u,H/u,R/1,R/0} for case b, and {C/u,H/u} for case c. 


lll. FAULT MODELS 


Most of the faults that are peculiar to the devices considered in this 
paper can be simulated using the PUF-PDF model of Figs. 7 or 9. 





Fig. 14—NMOs NOR gate. 
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3.1 PUF and PDF faults and enable line faults 


We shall consider first a special class of faults, where the PUF or PDF 
is enabled (disabled) when it is supposed to be disabled (enabled). This 
class is partitioned into four subclasses described in Table XII. Note 
that these are single faults. 


Table XII—Class of PUF/PDF faults 
Fault Sub- Fault Effect on Fault Effect on 


class PUF PDF 
I Disabling No effect 
I No effect Disabling 
I Enabling No effect 
IV No effect Enabling 


The fault effects on the output are described in Table XIII. The 
values marked with * are technology dependent and possibly unknown. 
Subclasses I and II cause a regular bistate gate to display tristate 
properties and a tristate device to be disabled when it is supposed to 
be enabled. Subclasses HI and IV may cause a conflict (overlap) under 
the appropriate input values. In some sense, this fault class blurs the 
difference between a regular bistate gate and a tristate device. For this 
reason, we can use the same set of logic values for faulty and fault-free 
circuit modeling, namely a TTL set (set 1 or 2) or a CMOS set (set 3), for 
both tristate and bistate devices. The only difference between the two 
cases is that the high-impedance state will not be produced during the | 
normal operation of a bistate device. 


Table XIII—PUF, PDF, and output values 
Fault Pull Up = Other Func- 


or Pull Down tion Output Value 
Yu =0 Yo=0 Y=H/* 
Yvu=0 Yo=1 Y=R/0 
Yu=1 Yo=0 Y=R/1 
Yu=1 Yp= Y=C/* 
Yp=0 Yu=0 Y=H/* 
Yp=0 Yu= Y=R/1 
Yp=1 Yu=0 Y=R/0 
Yo=1 Yu=1 Y=C/* 


Practically, the faults in each subclass can be obtained by shorted or 
open junctions in the PuF or PDF. We shall consider the example of 
Fig. 8 and summarize these subclasses of faults in Table XIV. In the 
context of concurrent fault simulation, the fault-injection mechanism 
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is extremely simple: a fault will be simulated if its effect on the gate is 
different from the fault-free circuit behavior, and this difference is 
measured over the set of possible logic values. 

Another class of faults can appear in tristate devices and concerns 
the enable line. An enable line stuck at 0 (stuck at 1) will cause the 
device to be permanently enabled (disabled). A permanently disabled 
gate can be modeled by an output “stuck at Z.” In the latter case, such 
a fault is characterized by a Z appearing at an output instead of a 
known logic value. The faults “enable line stuck at 1” and “output line 
stuck at Z” are equivalent. 


Table XI1V—tTypical faults in a CMOS NAND gate 


Fault Fault- 
Subclass Fault Example Input Values Free Fault 
I P, open m=0x%=1 1 Z* 
II N2 open m=lxm=1 0 Z* 
Ill P2 shorted xy =lx=1 0 a 
IV N, shorted 1 =0x%=1 1 a 


* With memory of previous logic value. 


3.2 PMOS and NMOS devices 


The structure of an NMOS (PMOS) device is similar to the cmos 
structure in that there is a pull-down function, but the pull-up function 
is a degenerate case of the cmos pull-up function, namely it is perma- 
nently enabled and serves as a resistor (Fig. 14). The function c = 
a + bis implemented in Fig. 14. 

We shall consider several fault modes, namely N, open, N; shorted, 
Ns open, and N3 shorted. The behavior of these four fault modes is 
represented in Table XV. The previous and present values of c are 
denoted by c(—) and c, respectively. The behavior of faults N; open, 
N, shorted, and N; shorted: is independent of c(—). However, if N3 is 
open, it is impossible to set c to a one, whereas the combination: 


c(—) =0 
a =0 
b =0 


produces a disabled output with new logic value equal to previous logic 
value (c = ZO). In spite of the behavioral differences, it is possible to 
model pMos and NMOS devices using the same set of logic values as for 
CMOS devices. 
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3.3 Input-open faults in CMOS gates 


It was shown earlier that a disabled cmos tristate device displays 
certain memory properties that could be modeled using additional 
logic values ZO and Z1. When an input to a CMOS gate is open, it is 
possible to produce these logic values in the faulty circuit. One method 
of modeling this is by setting the signal value at the site of the fault to 
a special value “propagating Zi’ (i = 0, 1, u), and propagating the effect 
to the gate output. Denoting the propagating ZO and Z1 by PZO and 
PZ1, respectively, we have the following conditions for the generation 
of these logic values at the site of the input-open fault: when the input 
changes from 1 or Z1 (0 or ZO) to 0(1), the faulty value of the input 


Table XV—Typical faults in an NMOS NOR gate 
















Open N, Short M Open N3 Short N3 
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Impossible 
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becomes PZ1 (PZ0). With the introduction of these additional logic 
values for modeling the fault, we need a method of propagating them 
through gates. Table XVI shows the NAND function whose inputs are 
from the set {0, 1, Z0, Z1, PZ0, PZ1}. Since we are considering single 
faults, four-entries in Table XVI are undefined. 

This modeling may be applied to the NAND gate of Fig. 8. Consider 
the fault, junction Ni open, and x; passing from 0 to 1 while x2 = 1. 
The fault-free output will pass from 1 to 0 and the faulty output will 
remain at the value Z1, meaning that the fault may be detected after 
the change of x, to 1, if it is possible to register Z1 as the value 1. 

Input-open faults in cmos gates can also be treated as special types 
of faults that may produce ZO or Z1 on gates outputs depending on 
gate type and present and previous input values. However, this would 
require treating input open faults on different types of gates differently. 
The proposed method presents a uniform way of inserting the effect of 
the input-open fault and only needs additional logic values during gate 
evaluation. These logic values introduced for modeling the fault do not 
themselves reach the gate output. 
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Table XVI—NAND function with propagating high- 
impedance states 


0 1 ZO Z1 PZO PZ1 
0 1 1 1 1 1 1 
1 1 0 1 0 Z1 ZO 
ZO 1 1 1 1 1 1 
Zl 1 0 1 0 Z1 ZO 
PZO 1 Z1 1 Z1 — _— 
PZ1 1 ZO 1 ZO — — 


IV. CONCLUSIONS 


A general model consisting of a pull-up function, a pull-down func- 
tion and an integrator is proposed for modeling TTL totempole, cmos, 
PMOS, and NMOS devices. It is shown that an accurate representation 
of the state of tristate devices and also certain bistate devices require 
not only the logic values but also impedance values. A set of 12 
combinations of impedances and logic values is proposed, each of 
which can be represented by a single value or by an impedance/logic 
value pair. Speed-storage trade-offs will determine the choice of rep- 
resentation. The set of logic values needed is shown to be a technology- 
dependent subset of the 12 combinations represented. The proposed 
model covers all the known tristate fault-free and faulty effects and 
does not require any additional modeling gates. 


REFERENCES 


1. R. L. Wadsack, “Fault Modeling and Logic Simulation of cmos and Mos Integrated 
Circuits,” B.S.T.J., 57, No. 5 (May-June 1978) pp. 1449-74. 

. R. H. Krambeck, private communication. 

. J. S, Jephson, R. P. McQuarrie, and R. E. Vogelsberg, “A Three-Value Computer 
Design Verification System,” IBM System J., 8, No. 3 (1969), pp. 178-188. 

. E. B. Eichelberger, “Hazard Detection in Combinational and Sequential Switching 
Circuits, IBM J. Res. Develop., 9, No. 2 (March 1965), pp. 90-9. 

. E. G. Ulrich and T. Baker, “The Concurrent Simulation of Nearly Identical Digital 
Networks,” Computer, 7, No. 4 (April 1974), pp. 39-44. 


ao -, Wr 


LOGIC SIMULATION MODELS 1287 


Copyright © 1981 American Telephone and Telegraph Company 
THE BELL SYSTEM TECHNICAL JOURNAL 
Vol. 60,, No. 7, September 1981 
Printed in U.S.A. 


Special-Information-Tone Frequency Detection 


By A. FEUER 
(Manuscript received January 20, 1981) 


The ability to distinguish recorded announcements from other call 
types is an essential part of the mechanized service-evaluation proc- 
ess. To do this, all announcements will have a special-information- 
tone prefix. A frequency detector—which is based on the correlation 
functions of the received signal—would be used to decide which 
announcement was triggered by a specific call attempt. This paper 
evaluates the performance of the frequency detector in the presence 
of additive noise and frequency shift induced by the announcement 
machine. The theoretical results, based on a calibration frequency, 
are very encouraging. To verify that use of this frequency is feasible 
in practice, an algorithm is proposed and its performance evaluated 
to show that it compares favorably to the theoretical one. 


1. INTRODUCTION 


As part of the process of evaluating the end-to-end performance in 
the telephone network, a sample of call attempts is evaluated and the 
attempts are classified into several categories, such as completed, busy, 
recorded announcements, etc. To mechanize this classification process, 
a machine must have the ability to distinguish between completed 
calls and recorded announcements. Current planning for the mecha- 
nized systems envisions the use of special-information-tone (SIT) pre- 
fixes which are to be attached to recorded announcements and can 
then be automatically recognized by the mechanized classifier (as well 
as alert the customer to the fact that he is listening to a recorded 
announcement). By choosing four distinct sITs, each representing a 
certain category of recorded announcements, the classifier will have 
the ability to distinguish between these categories as well as to recog- 
nize a recorded announcement in general. 

The siT is defined as a sequence of three consecutive tones. To get 
four distinct siTs five frequencies were chosen: fi<fp<fi<fi< fs, and 
each SIT consists of fi, fj, fs. The third tone has a fixed frequency and 
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will be used for calibration as is described later. The actual values of 
these frequencies are 904.5 Hz, 985.4 Hz, 1356.8 Hz, 1440.2 Hz, and 
1758.5 Hz. (The choice of these frequencies is mainly the result of 
constraints imposed by the cciTT definition of special information 
tones and tone generation and detection considerations.) To recognize 
which of the possible sits was received, the machine has to detect in 
the first tone whether it was fi or f2 and in the second tone whether it 
was f; or fz. This means that the classification process of the announce- 
ment categories is reduced to the two frequency-detection processes 
mentioned above. However, the planned direct recording of the siIT 
followed by the recorded announcement on the various announcement 
machines in the telephone network may introduce significant degra- 
dation into the reproduced sIT. In addition to additive noise, which is 
common to all signals in the network, frequency flutter and frequency 
shift of considerable effect on the reproduced tones may occur. The 
flutter effect, having an oscillatory nature, can be minimized by aver- 
aging properly the received data. In this paper, we report on our 
investigation of the combined effects of additive noise and frequency 
shift on the detection process. The detection scheme considered here 
involves use of correlation functions, and we have concluded that a 
reliable detection of sit frequencies is possible provided that a certain 
level of signal-to-noise ratio is ensured and the available data is 
properly used. It should be pointed out that in order to carry out the 
analysis we assume that the additive noise is white. We recognize, 
however, that in reality a colored additive noise may have a dominating 
effect on the error bounds. 

The structure of the paper is as follows. Section II presents the basic 
ideas of frequency detection via the correlation functions. In Section 
III, a white additive noise is considered to be present and the perform- 
ance of the detector in this case as a function of signal-to-noise ratio 
is evaluated. The effects of frequency shift are introduced in Section 
IV and the use of the calibration frequency to eliminate these effects 
is considered. The performance of the frequency detector as a function 
of the signal-to-noise ratio is evaluated and the results of using the 
calibration frequency are compared to the results when it is not used; 
a significant improvement is observed. 

The performance evaluations in Sections III and IV are based on 
using the likelihood-ratio test. This, because of the complexity of the 
expressions involved, is not practically feasible; however, these evalu- 
ations provide bounds on the performance of other schemes. In Section 
V, an algorithm for a detection process is presented. This algorithm is 
simple enough to be practical and its performance compares favorably 
to the theoretical one. 
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ll. CORRELATION DETECTOR 


The use of correlation functions of a received signal for detection 
purposes was devised as part of the overall classification process by J. 
E. Walls of Bell Laboratories.’ Some of the principles of the classifi- 
cation process that are relevant to our analysis are briefly described 
here. 

Let r(t) be a received signal observed for a period of length T. The 
correlation functions are defined by 


1 T/2 
C, == r(t)r(t + nd)dt n=0,1,2---, (1) 
T T/2 


where A < T (practically we want A to be small relative to T). 

We note that C> is in fact the signal’s average power in this time 
interval. 

Since we are interested in frequency detection, let us assume for a 
moment that the received signal is a pure frequency, namely 


r(t) = sin 2nft 


and derive expressions for the correlation functions in this case. Then, 


1 T/2 
C, == sin 2z7ft-sin 2nf (t + nA)dt, 
£ T/2 


or carrying out the simple integration results in 
c, = £98 2nfnd|_ sin 4rfT 

7 2 AnfT | 
The power term Co and the next two terms C;, C2, are used for the 
detection process. Since 


1 sin 4xft 
SS eee 
oe Anft 


C, = Co cos 2afA 
C2 = Co cos 2(27fA), 
if we normalize C, and C2 by the power Cy we get the following 


relationship 
Co Ci] 
[et] =f] — r 


Expression (3) means that all pure frequencies correspond to a point 
on a parabola in the (Ci/Co, C2/Co) plane. Clearly this parabola exists 
only for —1 = C,/Cy = 1 and —1 S C:/Cy = 1 (see Fig. 1) and Ci/Co, 


(2) 
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_ f= 1758.5 


f= 1356.8 Hz—~ 


~— f = 985.4 Hz 


Fig. 1—Correspondence between various frequencies and the parabola on the (C,/Co, 
C2/ Co) plane. 


C./Co are periodic functions of the frequency, so more than one 
frequency may correspond to the same point on the parabola. This 
means that for the detection process one must be concerned with 
the uniqueness of the points corresponding to the frequencies in- 
volved. For instance, if the two frequencies involved are f; = % A and 
fz = % A, both correspond to the point (0, —1) in the (Ci/Co, C2/Co) 
plane. This means that computing Co, Ci, C2 in this case is not sufficient 
to make the detection between these two frequencies possible even in 
the deterministic case (absence of noise). 

However, the choice of the parameter 1/A for a given set of frequen- 
cies involved in the detection process can ensure the necessary unique- 
ness. One way of doing it is to choose 1/A to be larger than twice the 
largest frequency to be detected, and this will be sufficient, as can be 
seen in Fig. 1. 

In the presence of noise, as is shown later, the point will shift from 
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the parabola towards the origin, while the origin itself corresponds to 
white noise for which both C, and C2 are equal to zero. 

In practice, instead of integration, the received signal is sampled and 
the sampled values are then used to get a good approximation of the 
correlation functions. Let r(t) again denote the received signal, sam- 
pled with frequency F, = 1/A, and N the number of samples done in 
the observation interval T. Then eq. (1) will be replaced by 

N-n-1 


Ch a N x iT i+ns (4) 


r=Pr : 
C= F, | 


Again, for a pure frequency, after some manipulations it can be 
shown that . 


where 


Ch = LS sia f isi 2 Epi 
n= de mt sin 77 
can be written as 
sin N 2a L 
n= 5 cos 7 N nos "7 2 
F, 
1 , cos NaF sin n Qn L 
+ 5 N cos(N — 1) 22 | ———_——__——. (5) 


sin 27 Le 


s 


With the assumption that N is large compared to n = 0, 1, 2, we may 
write 


; f 
N 22 — 
C ae 2 f i (N — 1)2 ca iF (6) 
= COS ee tt 7 eek ue 5 
sin 27 


so that the relationship eq. (3) holds here as well. In the Appendix, the 
values of the correlation functions computed by eqs. (1) and (6) are 
given for the frequencies used for sir. The error values introduced by 


going from eqs. (5) to (6) are also given in the Appendix for the same 
frequencies. 


Suppose now that a frequency, one of possible two, is transmitted. 
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If on the receiving end we make sure that the sampling frequency is 
large enough compared to the two possible frequencies, we ensure the 
uniqueness of the correspondence between these frequencies and the 
parabola defined by eq. (3). With this assurance, we can sample the 
received signal, compute Co, Ci, C2 according to eq. (4), and then use 
the values Ci/Cy and C2/Cy to decide which frequency was initially 
sent. 


lil. ADDITIVE NOISE 


In the previous section, frequency detection using correlation func- 
tions was described. As long as the received signal contains nothing 
but the pure frequency originally sent, the problem is clearly deter- 
ministic. However, in practice there are various degradations that 
affect the received signal. In this section, we are interested in evalu- 
ating the performance of the correlation detector in the presence of 
additive white noise. 

Thus, let the received signal be 


r(t) =A sin 2zfrt + n(t) k = 1, 2, 
where A corresponds to the signal-to-noise ratio and n(t) is a normal- 
ized white Gaussian noise. The sampled values are 
r,;=A sin & i + ni, (7) 


where 


nore hima =| 04 = an, 


k=1, 2, 1=0,1,2,---,N-—1, 


and n; are independent identically distributed (11D) variables with zero 
mean and variance one. 
Substitution of eq. (7) into eq. (4) results in 


1 N-n-1 
C= a y A’sin 8, sin(i + n)6, 
i=0 


N-n-1 N-n-1 
+ ¥} Alnisn»sin 16 + nsinfit+tn)é]+ Y nies (8) 
1=0 i=0 

Since eq. (8) contains random variables, Co, Ci, C2, and hence, the 
ratios C\/Co and C2/C> are also random variables with certain density 
functions. Once these density functions are known the detection proc- 
ess becomes a standard problem described in detail in various text- 
books on detection theory (see for example Ref. 2). We use the 
likelihood ratio test. To compute the threshold, we assume that the 
frequencies to be detected have equal a priori probabilities. The costs 
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are assumed to be zero and one for correct detection and error, 
respectively. With the knowledge of the density functions and the 
threshold, we can evaluate the performance of the detection process 
by computing the error probability. 

The first step then, is to develop the expressions for the density 
functions of C\/Co and C2/Co. However, in view of the complexity of 
eq. (8), rather than attempt to develop an exact expression, we will 
make use of a version of the Central Limit Theorem (Theorem 4.2.5 in 
Ref. 3) and develop approximate expressions for these density func- 
tions. 

Let us now denote 


2 N-n-1 
je 


Si=— > sin 26, sin(t + n)h. (9a) 

No 

N-n-1 
y?,= - Y [nisnSin 16, + nsin(i + n)6z]. (9b) 

i=0 

1 N-n-1 

Yno = WN > NMNi+n. (9c) 
Then, 

Ch = Si + Yi + Yoo. (10) 


With some manipulation, we can get more convenient expressions 
for S* and Y%;: 








A’ N-n } sin N8@. 
se = 6 eos a k 
9 cos n@, NV vy cost 1) ani 


2 


di 
+N cos(N — 1)6, Sus 


sin 0, 





, A N-n-1 
Yn = 2 n°? nb, ) nisin i), 


Aft N-1 
= NV » masini + n)+ Y ngsin(i — nyt 
=0 N-n 


Now, since N >> n we can write 





A 1 sin N@ 
St =— cos nb, | 1 - = cos(N — : 
9 cos n@; Vy cost 1); nd 
and 
A N-1 
Yr =2—cosn& Y nsin ib,, 
N = 
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or as 





Ae 1 sin N@, 
| ee pees as 6 
So 5 E Vy cos 1)6; sin 0 (11a) 
A N-1 
Yi =2— ¥ nin i6,. (11b) 
N io 
We get 
S* = Skcos nOy (12a) 
Y*, = Yi.cos nb. (12b) 
Denoting 
Ci 
xk = Ci 
C3 
and using eqs. (10) and (12) we get 
Sé+ Y6é O, + Yi 
pe lee tt Sake (13a) 
So + Yo + Yoo 
Sk + Yi 26, + Y& 
xk - [ 0 01 |Cos k Y 2 (13b) 


S§+ Y3+ Yh 


We observe now that xf and x are functions of the random variables 
Y4,, Y4o, Yt, Y2, each one of which is a linear combination of 1D 
random variables. The commonly used version of the Central Limit 
Theorem can be applied to conclude that the above four variables 
have Gaussian limiting distributions. This in turn enables us to use 
Theorem 4.2.5 in Ref. 3 to find the distributions of xf and x3. 

The first step is to compute the means and variances of Yiu, 
Yé&, Yio, and Y:. Using the fact that the 7,’s are 11D Gaussian random 
variables with zero mean and variance one and using eqs. (9c) and 
(11b) we can readily verify that if 


Yai 

yz 

pe 02 
eae 

Y bo 

then 

0 
E(y*) = |} 
FTG 
0 
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and 


4 
N S§0 0 0 
2 
0 700 
coy) = | [Y*— B(yHIEY*— ECV] = 0 0 x0 (14) 
1 
0 005 


As in Ref. 3 denote vector b: 


b = E[Y*]= 


oor © 


From eq. (14), the entries of Y*, being Gaussian independent random 
_ variables, are jointly Gaussian, as is VN Y* — b) with zero mean and 
covariance matrix 





4 
= gi 
WN 00 00 
2 
0 —00 
ae 
1 
T=N 0 0 we ; 
1 
0 ses 
0 0 N 
Now, let 
k k x* 
w= | (15) 
xX2 
and 
ol ax? ; 
0 = Nave |” 
so by eq. (13) 
cos 6; cos 20; 
ee —Skcos 6, —Sécos 26; 
Jn[se+i1fz | Sot 0 
0 S§+1 


All this establishes the preconditions for Theorem 4.2.5 in Ref. 3, 
which then states that the vector [x*(Y*) — x*(b)] has a Gaussian 
limiting distribution with zero mean and covariance matrix 
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R=9¢T¢o. (16) 
Since N in our case is quite large, we may write 


1 


ky = 
Pxt(x") = pzy ,(| fe) ~ On| Ral? 


exp| - ; [x — x*(b)/ Re [x - ren (17) 


where 


R= 1 (S§)?[2 cos’, + 1] + 2Sé[2 cos?6, + 1] + 1, 
+ N(Sk + 1) 2S%cos 6,cos 20,(S% + 2), 
2Sécos Oz cos 26,(S% + 2), 
(S5)?[2 cos?20, + 1] + 2S%[2 cos?20, + 1] +1 
and 





E(x") = x*(b) = Sé cos 0; | 
Sk +1] cos 26, |’ 

From eqs. (lla) and (17), we can readily observe the way the signal- 
to-noise ratio, A, affects the density function. If A gets very large, Sé 
becomes very large; then the mean value of x* approaches the parabola 
and the entries of R, become very small—namely, we approach the 
deterministic case described in the previous section. On the other 
hand, if A approaches zero, S§ goes to zero as well and the mean value 
approaches the origin. This is understandable since as A gets smaller, 
the effect of the noise gets larger, dominating the signal; and, being a 
white noise, the cross correlation functions C; and C2 approach zero. 

In Figs. 2 and 3, we see the effect of A as described above for two 
pairs of frequencies that were selected for the sir. The sampling 
frequency is 4000 Hz. It should be noted that the mean values for each 
frequency are on a straight line as is obvious from eq. (17). In the 
figures, for every A, one equal-probability contour is drawn around the 
mean and the fact that these contours get smaller as A gets larger is 
expected since the entries of R; are getting smaller as was pointed out 
earlier. 

Once we have eq. (17) the detection is straightforward. With a priori 
probabilities and costs as described earlier, the measured data is used 
to compute the point in the (Ci/Co, C2/Co) plane and test which 
conditional density has higher value at this point. Then decide on the 
corresponding frequency as the one that was sent. More details of this 
procedure are described in Ref. 2. However, using the expressions we 
have for the density functions the process can be somewhat simplified. 
From eq. (17) and Ref. 2, 
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02 04 06/Cij/c, 





Fig. 2—Dependence of the density functions P(x| f) on signal-to-noise ratio for 904.5 
and 985.4 Hz (dashed lines are the equal likelihood points). 


fi 
Dx 7,(%| fi) = Px (| fr), 
2 


which means that if the inequality holds in one direction, we decide on 
hf, and if in the other, /2. This is equivalent to 


In| R2| — In| Ril + [x — x7(b) )'R2"[x — x7(b)] 
fi 
— [x — x°(b)/Ri"[x — x*(b)] = 0, (18) 


which with an equality sign is a line in the (Ci/Co, C2/Co) plane. The 
dotted lines in Figs. 2 and 3 correspond to the above-mentioned line. 
On one side of these lines one density function has higher values, on 
the other side the other density function. The detection then consists 
of deciding on which side of this line the measured point falls. 
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(c) 


Fig. 3—Dependence of the density functions p(x| f) on signal-to-noise ratio for 1356.8 
and 1440.2 Hz (dashed lines are the equal likelihood points). 


To evaluate the performance of this detector, we compute the error 
probability. Let B,, k = 1, 2 be the part of the plane in which 
Px\f,(x| fe) has the larger value. Then the probability of error will be 


Pr(e) =5 | 
B. 


Figure 4 shows the error probability as a function of A for, again, 
the two pairs of frequencies of interest in the sits 904.5 Hz, 985.4 Hz, 
and 1356.8 Hz, 1440.2 Hz, and sampling frequency of 4000 Hz. Since, 
in general, signals transferred by the telephone network result in 
signal-to-noise ratios higher than 8 dB, the results here are encouraging 
for both frequency pairs. It is interesting to note that the higher 


Dx f(x| fildx + | 


By 


Px| f,(x| Ade ~ (19) 


2 
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frequency pair results in a better performance. Observing Figs. 2 and 
3, one can see the reason for this difference in performance; the density 
functions for every signal-to-noise ratio are better separated in the 
higher pair. It turns out, that for every choice of a sampling frequency 
and difference in value, some pairs of frequencies—and not necessarily 
the lower valued frequencies—are more detectable than others. 


IV. ADDITIVE NOISE AND FREQUENCY SHIFT 


The use of announcement machines as tone generators will introduce 
two primary types of degradation, frequency flutter and frequency 
shift. In this section, we attempt to analyze the effects of the frequency 
shift with additive noise on the performance of the correlation detector. 
It is assumed that the flutter effect is eliminated by averaging a 
number of observations of each received tone. 

The complexity of the expressions involved in the analysis here 
makes closed-form results very difficult if possible at all. For this 
reason, digital computer calculations were used as the main tool in the 
analysis. 

Introducing the frequency shift, we get the expression for the re- 
ceived signal 

(1 + na) f 


r(t) = A sin Lage + n(t), 


904.5 Hz 


Pr (€) 


— ~ 985.4 Hz 





1356.8 Hz / 
1440.2 Hz 








A IN DECIBELS 


Fig. 4—Probability of error when no frequency shift is considered. 
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where na is assumed to be a random variable with extended beta 
distribution. This particular distribution is general enough to include 
many possibilities and agrees with the physical properties of na 
(namely —1 = ng <= 1). So na’s assumed density function is 


— {BL + na) "(1 —na)*" —s for nal <1 
Pn,{na) = 0 elsewhere’ 
where 
atB+1 
7 E T(a + B) 
2 T(a)P(B) 





and a,8 = 1 are the two distribution parameters (in our computations 
we chose a = B = 10 which fits reasonably the little data available for 
Na). ‘ 

The presence of this additional degradation causes the expressions 
developed in the previous section to be conditional on ng. This means 
that now, rather than having an expression for p,) ;, ("|"), we have an 
expression for pzjn,,/,(_ |") or using eq. (17) we may write 


1 
27|R.|'” 


exo| -3 [x — E{x"}]Ri'[x - ta) | (20) 


Px\ ngfplX | Na, fx) = 


where the expressions for R;, E{x*}, S§ are as before [see (1la) and 


(17)] and 6, = an LL Mad fe, 


Since for the likelihood ratio test p, ;,('|')is needed, we can proceed 
to compute it using the relation 
1 


Px, (x| fx) = | Px\ngf,A(X| May fe)Pn,(ta)dna. 
-1 
With this and eq. (19), the performance of the detector can be evalu- 
ated for this case where both additive noise and frequency shift are 
present. In Fig. 5 the probabilities of error in detection as a function of 
the signal-to-noise ratio are presented. Comparing these results to Fig. 
4 makes it clear that the performance of the detector deteriorates very 
significantly when frequency shift is present. Even for a very high 
signal-to-noise ratio the frequency shift induces considerable error. 
The effect of the frequency shift alone, which provides lower bounds 
on the error probabilities in Fig. 5, can be calculated as follows. When 
A gets very large the received signal becomes a pure tone with a shifted 
frequency. This, however, implies that in this case it will be sufficient 
to consider only x; = C\/C> for the detection. Since x; and ng are 
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__904.5 Hz 
— 985.4 Hz 


Pr(€) 


1356.8 Hz_ — 
1440.2 Hz 





—6 -—4 -2 0 2 4 6 8 10 12 14 16 
A IN DECIBELS 


Fig. 5—Probability of error when frequency shift is considered but no calibration 
frequency is used. 


related through 
fi 


x1 = cos 27(1 + na) F? 
the knowledge of pn,(') makes the calculation of p.,);,(°|°) for each 
frequency f;, straightforward. In Fig. 6 the density functions for the 
two pairs of frequencies are drawn and the thresholds for the detection 
in this case are the intersections of the density functions (also pointed 
out in the figure). 

The technique to overcome the frequency shift is based on using the 
third tone for calibration in the detection of the first two tones. This 
means that a fixed tone is sent, processed in the receiving end, and the 
knowledge of its exact value and the corresponding measurements can 
be used to improve the detection of the first two tones. 

To be more precise, let f; be the frequency of the third tone and 
x° =[Ci/C%, C3/C$] computed from the corresponding measured data. 
Then through eq. (20) we know that 


1 
Px3|\ nfl? | Na, fs) a 2Qa| Rs] 


exp \- 5 [x? — E{x3)R3 [x3 - e301 
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P (1 leg| fk) 





_-— fy = 905.4 Hz 






fo = 985.4 Hz ~ ~ 





-1 -0.5 0 0.5 1 
FOR THE LOWER INTERVAL FREQUENCIES 


(a) 


Pp (Cr/cq| fk) 





—1356.8 Hz 
1440 Hz—” cal 


| 
| 
| 
| 
| 


Ci/¢ 
-1 0.5 0 0.5 1 9 
FOR THE HIGHER INTERVAL FREQUENCIES 
(b) 


Fig. 6—Density functions required for detection when frequency shift is considered 
but the additive noise is ignored (A > ©). 


where the expressions for E{x*} and R3 are as in eqs. (11a) and (17) 
with f3 substituted for f,. Using this we can improve on our knowledge 
of the statistics of ng by computing pn,|x3,,(.|.,-) through the relation- 
ship 


P33 | nal (x? Nd, f3)Pn (na) 
Png\xiflna| x’, fal =z l"afs | d 


Pxd|nghlX* |Na, fa]Png(na)dna 
-1 
This improved data on ng can then be used to calculate px) x35, 


(. | seve) 
1 
Prizeg zl x|x'; fe] = | Px | ngflX | Na, fe)Png\x°,4[Na|x°, faldna, 


Sa) | 
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which in turn will be used for the one process, namely 
Px e41x|x", fr, A = = Pei frfLx|x*, fa, fe). (21) 


In Figs. 7 and 8, we present—for this improved detection process 
with 1758.5 Hz as the calibration frequency—the error probabilities as 
a function of the signal-to-noise ratio. Note that the sampling fre- 
quency 4000 Hz is larger than 2 X 1758.5. Comparing this to the results 
without the calibration frequency, we observe considerable improve- 
ment, whereas the results computed with no frequency shift present 
provide a lower bound on error probabilities (or an upper bound on 
the improved detector’s performance). 

The case when A — o (i.e., when the additive noise becomes 
negligible) is again of special interest but very simple. Since now with 
the knowledge of both x and fs the shift can be exactly computed, 

B 


x? = cos 2a(1 + na) 


or 


na= 5 cos” 1y? — de 


2afs 


— ao 
~ WITH SHIFT AND NO 
\\ CALIBRATION FREQUENCY 
| 


/ 


ve 


Pr(e) 


_WITH SHIFT AND 
~~ CALIBRATION FREQUENCY 





A IN DECIBELS 


Fig. 7—Error probabilities with frequency shift and calibration frequency—lower 
interval frequencies. 
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WITH SHIFT AND NO 
aa ila FREQUENCY 


7 
2 


Pr (€) 


_.WITH SHIFT AND 
CALIBRATION FREQUENCY 





-6 —4 -2 0 2 4 6 8 10 12 14 16 18 
A IN DECIBELS 


Fig. 8—Error probabilities with frequency shift and calibration frequency—higher 
interval frequencies. : 


and the detection of the first two tones becomes deterministic. The 
only source of problems in this case is the question of uniqueness of 
cos'x?. This can be overcome by considering both cos"’x} and 27 — 
cos ‘x? and with probability one, one of them will recover one of the 


two possible frequencies. 


V. A SUGGESTED PRACTICAL CALIBRATION ALGORITHM 


In the previous section, we described how the use of a calibration 
frequency can theoretically improve the detection of tones that are 
affected by frequency shift and additive noise. However, practical 
computation of the density functions required in eq. (21), or the 
separating line (where the two density functions are equal) is impos- 
sible. Hence, an algorithm is required that is both compatible to the 
theoretical exact density function separation and simple enough to be 
practically implemented. Such an algorithm is presented and its per- 
formance evaluated. 

Let us first review the available data to be used for the detection. 
The received sequence of three tones is sampled and for each tone the 
correlation functions are calculated repeatedly and averaged to elimi- 
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nate frequency flutter effect. Then the cross-correlation functions are 
normalized to give the three vectors 


1 Ci Ci 
Co C C3 
gles Wes = 
C2 C3 C3 
Co C C3 


corresponding to the three tones. Whereas the frequencies resulting in 
x! and x”, respectively, are not known (in each case they can be either 
one of two possible frequencies), the third one is known to result from 
1758.8 Hz. The idea is to use this knowledge to get an estimate of the 
frequency shift to help in the decision process of the first two tones. 
The suggested algorithm (see Fig. 9) is as follows: 

(i) Draw a line from the origin through x°; it will intersect the 

parabola at 


Cle. 






5(1356.8 + 1440.2)—— 


UPDATED__ ~ 
MIDFREQUENCY 2 


~__UPDATED 
MIDFREQUENCY 1 


7 
1 2 
7 (904.5 + 985.4) — 


Fig. 9—Geometric interpretation of suggested algorithm. 
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2 
x) 
8 
Mea} noce 


where 
x3 + v[x3]? + 8[ x2]? 
2x3 ; 


(ti) Use <3 to compute the values 


meee) \ = 1,2, (23) 


(22) 


where 


pez 904.5 + 985.4 
2% 1758.5 


= 0.5374 cos 1x3 (24a) 


-123 
x1 


and 


4, — 1356.8 + 1440.2 | 125 
9X 1758.5 : 


= 0.7953 cos £3. (24b) 


(iit) Draw the lines from the origin to the points (ai, a): If x’ is 
counterclockwise away from the corresponding line, the ith tone is of 
the lower frequency, and if it is clockwise away, it is uf the higher 
frequency. 

The motivation behind this algorithm is quite simple. We have 
observed earlier that the additive noise effect is to shift the mean value 
of the pair (Ci/Co, C2/Co) towards the origin, whereas the frequency 
shift causes this pair to move along the parabola. The first step in the 
algorithm can be viewed then as isolation of the effect of the frequency 
shift. The point £* on the parabola is regarded as the result of the 
original frequency 1758.5 Hz, together with some shift that can now be 
estimated using the relationship 


1758.5 
~3 2 A 
Xi = cos 29 00 (1 + 7a) 
or 
P 4000 psi 
OO — —_ 1. 
oe see 
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This estimate is then used to update the midfrequencies for each tone, 
%(904.5 + 985.4)(1 + nag) for the first and %(1356.8 + 1440.2)(1 + nga) 
for the second. The lines that connect the origin with the points 


m(904.5 + 985.4) 


oe (Le 
cos aod (1 + 7a) 
277(904.5 + 985.4) z 
aa 
and 
(1356.8 + 1440.2) . " 
cos ang, (1+ Na) 
27(1356.8 + 1440.2) (1 + ny) 
4000 oe 


provide the threshold lines for the detection of the frequencies of the 
respective tones. In Fig. 9 the algorithm is described geometrically. 

In Figs. 10 and 11 the performance of a detector using this algorithm 
is compared to the performance when exact separation is assumed. It 
is quite clear that the use of the algorithm results in a performance 
that is very close to the theoretical one. 


SUGGESTED 


EXACT ALGORITHM 
SEPARATION } 
\ 


~ ae 


Pr(€) 











A IN DECIBELS 


Fig. 10—Comparison of error probabilities for suggested detection algorithm to the 
exact separation—lower interval frequencies. 
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ALGORITHM 
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Fig. 11—Comparison of error probabilities for suggested detection algorithm to the 
exact separation—higher interval frequencies. 


There are two difficulties in the described algorithm that will affect 
its performance. The first is because positive shifts above approxi- 
mately 14 percent will result in frequencies higher than the critical 
one—F’,/2 (see Appendix)—for the third tone, and the points on the 
parabola are no longer uniquely related to their corresponding fre- 
quencies. This difficulty is inherent to the correlation function ap- 
proach and can be eliminated only if a significantly higher sampling 
frequency is chosen (approximately 7000 Hz). However, if we assume 
that the shift is always less than the critical one, even for the proposed 
sampling frequency, the effect does not seem to be significant since in 
the model we have chosen for the frequency shift the probability of 
having shifts higher than 14 percent is quite low. The second difficulty, 
which is inherent in the proposed algorithm, arises when the resultant 
X? is less than —1. In this case, we propose simply to take it equal to 
—1 and again argue that the probability of this happening is very low 
even for small signal-to-noise ratios. Altogether, both difficulties, if 
treated as is suggested above, do not seem to affect the performance 
of the detection algorithm. 


VI. CONCLUSION 


In this paper, we addressed some of the problems in detecting 
recorded announcements encoded via the special information tones 
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(sit). In particular, we discussed problems that arise when additive 
noise is present. We have assumed that the frequency flutter effects 
are eliminated by averaging several observations, and investigated in 
detail only the additive noise and frequency shift effects. 

Our results support the conclusion that by properly using the infor- 
mation on the frequency shift, its effects can be made almost negligible 
and under these conditions high-performance SIT detection can be 
achieved. 

The performance evaluations presented here make explicit use of a 
certain assumed model for the noise and frequency shift; however, the 
detection algorithm, which is proposed in Section V, is independent of 
any such assumptions. The performance of this algorithm is compa- 
rable to that theoretically achievable, and thus this algorithm is 
proposed for implementation in the siIT classification process. 
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APPENDIX 


Approximations Used in Computing the Correlation Functions by Sampled 
Data. 


f 0 
(Hz) (rad) CC CC, CCz Co C; C2 Ei E2 


904.5 1.42 0.5 0.0747 —0.4778 0.4971 0.0743 —0.4749 0.0001 0.00004 
985.4 1.55 0.499 0.0114 —0.4985 0.4982 0.0114 —0.4976 0.0008 0.00003 
1356.8 2.13 0.4993 —0.2654 —0.2171 0.499 -—0.2653 —0.217 0.0008 —0.0008 
1440.2 2.26 0.4993 —0.3184 -—0.0932 0.4997 -—0.3186 -—0.0933 0.0003 -0.0004 
1758.5 2.76 0.5004 —0.4649 —0.3632 0.496 -—0.4608 0.36 -—0.0019 0.0034 


Continuously computed correlation functions [see eq. (2)]: 


CC, = ; eosGiby E ann | 





Anft 
Approximated correlation functions [see eq. (6)]: 
sin 9 


at eee = sin (N@) 
C, = 5 cos(n6) {1 v cos[(N — 1)6] | 


Differences between discretely computed correlation functions [see eq. (5)] 
and their approximated values. 


sin(n@) 
sin 6 





1 
E,= oN cos[(N — 1)@] cos(N@). 


[Note: F, = 4000 Hz; N = 167; T = 4% seconds and @ = 2m (f/F;)]. 
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Annoying echoes generated by impedance mismatches can occur in 
long-distance telephone networks. An echo canceler controls echo by 
synthesizing an echo replica and subtracting it from the actual echo 
on the return path. For the echo canceler to work properly, it must be 
able to distinguish between desired signals and echoes. One instance 
where this can be difficult is during the common-channel interoffice 
signaling (ccIs) continuity check between four-wire switching offices. 
On ccis trunks the voice and the signaling are routed separately. To 
ensure a satisfactory transmission path, a voice path continuity check 
is conducted before call setup. The check is performed on a loop basis 
by sending a check tone and looping it back at the distant office. 
Since the echo canceler adaptive filter memorizes information from 
the last previous call, it may partially cancel the looped-back tone. In 
this paper we study the effect of the echo canceler adaptive memory 
on the ccis continuity check. The analytical and experimental results 
indicate that the occurrence of continuity check failure caused by the 
presence of an active echo canceler is so infrequent and insignificant 
compared to the existing statistics of all other pertinent failure 
mechanisms. Thus, disabling the echo canceler during the ccIs con- 
tinulty check is unnecessary. 


|. INTRODUCTION 


An echo may be produced in a transmission system whenever there 
is an impedance mismatch. This impedance discontinuity can cause a 
significant portion of the transmitted signal energy to be reflected 
toward the signal source over an echo path. Noticeable echoes consti- 
tute one of the most serious forms of impairment in telephone channels. 
Its subjective annoyance increases with both echo amplitude and 
propagation delay. The increasing use of satellite circuits for domestic 
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and international calls has made the control of echoes even more 
important. The traditional approach of echo control using via net loss 
or echo suppressors is not acceptable for full-hop satellite circuits, i.e., 
circuits carrying both directions of transmission via satellite.’” A better 
approach is to use echo cancelers that control echo by synthesizing a 
replica of the echo and subtracting it from the returned signal.” The 
realization of the single-chip integrated echo canceler has made its 
widespread deployment appealing.* 

Figure 1 models the essential elements of the echo path to show how 
echo cancelers work. When there is far-end speech x(¢) but no near- 
end speech vu(f), the internal registers of the echo canceler adaptively 
update the estimate of the echo path impulse response A(t) to form an 
echo replica y(t) that is subtracted from the real echo y(t). When the 
near-end speech is detected in the presence of the far-end speech 
(double talk) the speech detector inhibits further updating, but the 
echo canceler still tries to cancel the echo contained in y(t) by using 
the most recent estimate of the echo path impulse response. This 
property of continued echo canceling during double talk is a nice 
feature of the echo canceler. However, its effect on the common- 
channel interoffice signaling (ccIs) continuity check between four-wire 
switching offices has caused some concern.” 

Common-channel interoffice signaling is a system for exchanging 
information between processor-equipped switching systems over a 
network of signaling links. All signaling data for call setup and take- 
down, as well as network management signals, are exchanged by these 
systems over the signaling links instead of being sent over the voice 
path. Thus, a continuity check for voice path assurance (VPA) is 
conducted whenever a ccis trunk is selected to switch a call forward. 
This check not only ensures a satisfactory transmission path but also 
precludes billing for an otherwise undetectable faulty connection. 
Between four-wire switching offices the vpA check is performed on a 
loop basis by connecting a single-frequency transceiver at the origi- 
nating office and looping the transmission pairs at the distant office. 
The check is considered successful when the vPa tone received at the 
originating office is within acceptable transmission and timing limits. 

Figure 2 shows the vpa check for the satellite circuit. At the distant 
office the same notation x(¢) is used for the received and the transmit- 
ted signals because of the looping configuration. The echo canceler at 
the distant office will detect double talk during the vpa check. If the 
echo canceler is not disabled, it may partially cancel the check tone 
according to the last estimate of the echo path impulse response 
because of its continued canceling property during double talk. There- 
fore, the introduction of the echo canceler into the long-distance 
telephone network without proper disabling control will have some 
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Fig. 1—Use of echo canceler to control echo. 
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Fig. 2—Echo canceler and ccis continuity check. 


influence on the ccIs continuity check. Voice path assurance failures 
will lower the performance indices of the affected offices and generate 
additional maintenance activities. Unfortunately, since call processing 
information is not readily available, no satisfactory solution has been 
found such that the echo canceler can be disabled at the desired 
moment. The seriousness of the echo canceler effect on the vpa check 
is studied below by evaluating the probabilities of various possible 
events caused by vpaA failures due to the use of the echo canceler. 
Section II describes the circuits and facilities in the VPA path and how 
they are modeled mathematically in the derivations given in the 
Appendices. Section III presents analytic and experimental results and 
discusses the impact of echo canceler on the cCcIs continuity check. 


ll. MODELING OF VPA CHECK 


In this section, the discussion and the derivation will be restricted to 
No. 4 Ess offices and satellite circuits because the results can be 
similarly derived for other facilities. 

The vpa check for ccis-equipped No. 4 Ess offices is a two-step test 
in a specified time interval. The first step is to detect at the proper 
receive level the transmitted vPpaA tone which is looped back at the 
distant office. The second step is to stop transmitting the vPA tone 
and detect its disappearance at the receiver. At present the specified 
time interval for completing the test is set anywhere between 2 to 3 s, 
depending on the specific office. At the beginning of a vPpaA check 
during call setup, the originating processor of the switching machine 
starts the timer and attaches a 2010-Hz transceiver to the selected 
satellite trunk concurrent with sending an address message through 
terrestrial ccis links identifying the trunk to be looped back. The 
distant office, upon receipt of the message, connects the receive side of 
the trunk to the transmit side through a zero-loss loop. The originating 
office checks the level of the returning tone to verify that transmission 
loss is within acceptable limits. If the returning tone has acceptable 
level, the originating office stops sending the vPA tone and verifies that 
the receiver measures below a release level. A vPA failure is generated 
if for any reason the above two steps are not completed before the 
timer times out. 

The specification of the transceiver levels is given in Table I. 
Frequent and periodic tests are performed in switching offices to make 
sure that the transceiver levels are within specification. Table I also 
gives the test requirements which are stricter than the specification.°® 
Since the test only gives pass or fail indication, it is assumed that the 
various levels tested are uniformly distributed within the test limits. 
Thus, relative to each other, the transmit and the detect levels of the 
VPA tone are assumed to be uniformly distributed in (—1, 1) dBm and 
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Table |—Common-channel interoffice signaling continuity check 


. levels’ 

Test Specification Test Requirement 
Transmit —12 + 1dBm0 ~—12 + dBm0 
Accept ~18+n=NdBm -18.6+n=NdBm 
Fail Ns -22+ndBm N <= 21.4+ndBm 
Release N = -—-27+n dBm — 


1 Where N is the absolute power level of the va tone and n is the relative power at 
the transceiver with respect to the zero transmission level point. 


(—9.4, —6.6) dBm, respectively. The release level in Table I is not 
modeled since a properly working echo canceler should not affect it. 

In addition to the transceiver level variations, the vPA tone level is 
affected by the two-way satellite trunk loss, the check loop, and the 
echo canceler at the distant office. The satellite trunk is a class of 
intertoll trunks with 0 dB inserted connected loss which is specified to 
be normally distributed with a standard deviation of around 0.7 to 0.8 
dB. The check loop is specified to have a loss of 0 + 0.1 dB. This 
variation is small and will be ignored in the following study. The effect 
of the echo canceler requires detailed consideration. 

In discrete-time notation, the near-end speech detector inhibits 
updating the echo canceler registers at time & if 


yk-D)>% max x(k-l—n) for 0<1< 255, (1) 


where the notations are given in Fig. 1.°> At an 8-kHz sampling rate, 
0 <n = 127 implies using a 16-ms tapped delay line adaptive digital 
filter to model the echo path, while 0 = 7 = 255 means that once near- 
end speech is detected, the detector continues declaring its presence 
for the hangover time of 32 ms. The factor of one half is based on the 
assumption that there will be at least a 6-dB loss through the hybrid. 
During the continuity check, the vPA tone is looped back so the near- 
end speech detector finds that eq. (1) is satisfied and thus inhibits 
updating the registers. The echo canceler will partially cancel the vPA 
tone according to the register settings determined by the echo path 
return loss (EPRL) of the last call. Thus, the degree of cancellation is 
a random variable which is a function of the EPRL. The EPRL distri- 
bution of the whole telephone network is not known. The only available 
EPRL data, as shown in Fig. 3, were measured at the Pittsburgh 
Regional Center in the second half of 1976 during the field evaluation 
of domestic satellite.’ These measurements are not inconsistent with 
previously reported echo path and intertoll trunk loss.’ Although the 
loss distribution was derived from the estimated echo path frequency 
response as averages over 500 to 2500 Hz, it will be assumed to be 
the loss distribution for the vPA tone. The density function shown in 
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Fig. 3—Distribution of near-end echo path loss, 500 to 2500 Hz flat-weighted. 


Fig. 3 is consistent with a normal distribution, except the small 
secondary peak at about 5 dB. The following derivation turns out to 
be independent of the loss distribution below 6 dB. 

The above modeling of the vpA path can be used to predict VPA 
failures which may generate various kinds of maintenance activities 
affecting trunks and carrier groups. There are 12 trunks in each carrier 
group. If a trunk in a carrier group fails a vpA check, the call is 
reattempted on another trunk in a different carrier group, while the 
first carrier group is temporarily taken out-of-service. Several seconds 
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later, a second trunk in the locked-out carrier group is randomly 
selected and checked for continuity. If the second trunk also fails, a 
software carrier group alarm (SCGA) is generated for maintenance 
actions. If the second trunk passes the check, the temporary lockout 
on the first carrier group is released and the initially failed trunk is 
rechecked. If it fails again, a single trunk lockout (STL) is generated for 
maintenance. The probability of getting a vPaA failure is derived in 
Appendix A. The probability of generating an STL is given in Appendix 
B. The probability of generating scGA cannot be analyzed since the 
two failed trunks which share identical facilities from the switching 
offices to the satellite earth stations may be dependent. Thus, only the 
upper and the lower bounds of the probability of scGA are obtained in 
Appendix C. 


ill. DISCUSSION OF RESULTS 


Figures 4 to 7 show the results obtained in the appendices. During 
the vpa check, the estimated echo <(¢) in Fig. 2 is also a single- 
frequency tone at 2010 Hz. Its amplitude and phase depend on the 
EPRL of the last call. For the worst-case, 6-dB EPRL, its effect on the 
looped back tone x(t) varies from 3.5 dB (if the tones add construc- 
tively) to —6 dB (if the tones add destructively). With the EPRL 
distribution given in Fig. 3, the effect of the echo canceler on any VPA 
tone level follows the density function shown in Fig. 4. This density 
has a mean value slightly larger than zero; that is, the expected effect 
of the echo canceler is to increase the vPa tone level slightly. This is 
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Fig. 4—Probability density function of the effect of the echo canceler on the vpa tone 
level. 
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Fig. 5—Probability of vpa failure. 


to be expected since the resultant VPA tone contains the powers of two 
sinusoids x(t) and <(t) with random relative phases as opposed to the 
single sinusoid x(t) without the echo canceler. However, the reducing 
effect on the vPA tone level has a greater impact because it can be as 
large as —6 dB. This is evident in Fig. 5, which shows that using echo 
canceler appreciably increases the probability of vPA failure. However, 
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Fig. 6—Probability of single-trunk blockout. 


the calculations assume that all circuits are working properly and 
without outside interferences. Voice path assurance failures can be 
generated by carrier glitch, fading, and many other causes. Limited 
sampling performed during the spring of 1980 in several No. 4 Ess 
offices not equipped with an echo canceler yielded vpPaA failure rates 
ranging from 0.5 < 107° to 2.7 x 10°*, with a sample mean of 5.6 X 
10°°.° These are of the same order of magnitude as the calculated 
probability of vpa failure with echo canceler, assuming the satellite 
trunk standard deviation is between 0.7 to 0.8 dB. Therefore, if the 
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measured vpa failure rate for a group of trunks is 2 x 10°, it will be 
probably be around 4 X 10° after echo cancelers are installed on these 
trunks. The calculated probabilities of stB and SCGA as given in Figs. 
6 and 7, respectively, also do not appear to be excessive, although there 
is no field statistics for comparison. 

An experiment was performed on the full-hop satellite circuits 
between Atlanta, Georgia, and Cedar Knolls, New Jersey to see the 
effect of not disabling the echo canceler prior to the vPA check. The 
data indicated that there was no vPA failure for well over 100,000 calls. 
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This result is better than the calculated probabilities and it can be 
easily explained. The calculated numbers are based on the specifica- 
tions of the trunks and circuits. The standard deviation of the satellite 
trunks appears to be smaller than that of the intertoll trunks. Fur- 
thermore, limited measurements of the actual transceiver levels indi- 
cated that the transmit and the detect levels can be modeled as 
uniformly distributed in (—0.25, 0.25) dBm and (—8.75, —7.75) dBm, 
respectively.’° All these factors help the performance of the vpa check. 

In conclusion, the introduction of echo canceler to the satellite 
circuits may generate additional vPaA failures. However, in present long 
distance networks, they are not significant compared to existing VPA 
failure statistics of the field. Thus, it does not appear necessary to 
disable the echo canceler prior to each CCIS continuity check. 


APPENDIX A 


The probability of vpa failure is derived in this appendix. Let / be 
the EPRL shown in Fig. 3; that is, / is normal with mean p = 15.37 and 
standard deviation o = 4.4, except for 7 < 6 dB. By definition, for a 
single frequency tone, 
power of x(z) 


ae power of y(t) 
amplitude of x(¢ 
= 2 ee 
WAGE amplitude of y(¢ 
amplitude of x(¢ 


= 20 log ——______ 
108 amplitude of y(t 


? 


where the last step follows for / > 6 dB. For the loopback configuration 
in Fig. 2 during the vpa check, then 


1 = 20 log =, 1>6 dB, (1) 


where 
_ amplitude of x(t) 
~~ amplitude of x(t)’ 


The derivation below assumes that eq. (1) is an exact equality. Since 
l> 6 dB, b is distributed in (0, %). Let J;,i = 1, 2, --- , denote the EPRL 
of the ith previous call using a specific echo canceler. The /;’s are 
identically distributed as /, and for most practical situations, independ- 
ent random variables. Clearly 6 is a function of the /,’s. For instance, 
if 1, > 6 dB, b is determined solely by i. If 1: < 6 dB, the near-end 
speech detector would declare double talk and no register update 
would take place during the last call. Then b would be independent of 
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l, and become a function of only J;, i = 2, 3, --- . Thus, the distribution 
function of b is given by 


P{b <c} = P{107™ < c} + P{l: < 6}P{10-”” < c} 
+ P{l, < 6}P{k < 6}P{10-"/" <e} +--- 
=(1+(1—g) + (1—g)’ + ---)P{10-"™ < c} 





1 1 
= Pt > 20 log | 
g c 
es 2 gt wtgy, (2) 
& 20log 1/e V200 
where 
g = P{l>6} 
ee | 
= exp[—(J — p)?/207] dl. 
| V270 
Differentiating eq. (2) gives the density function of b as 


20 1 1 
(b) =~ exp[—(20 log b + p)?/207], 0<b<-. (3) 
fe V2rog In 10 8 2 


The transmitted vPA tone a(t) can be written as 
a(t) = V2A sin wt. 


The transmit level in dBm 
A? 
x1 = 10 log tos 

is uniformly distributed in (—1, 1), as described in Section II. Since the 
satellite trunk loss in dB is normally distributed, it can be considered 
as gain for convenience. Assume the satellite trunk gains in the 
transmit and the receive directions are denoted by gi and ge, respec- 
tively. Let 


x2 = 10 log gi 
and 
x3 = 10 log gi. 


The terms x2 and x3 are independent normal random variables with 
zero mean and standard deviation p. Thus, using the notation of Fig. 
2, 
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x(t) = V2Ag; sin wt, 
e(t) = V2Ag, sin wt — V2Agib sin(wt + 6) 
= J2Ag,v1 + b? — 2b cos ¢$ cos E + tan™ (Ga *) | 
b sin } 
where ¢, the relative phase between the loopback signal x(t) and the 


echo estimate X(t), is assumed to be uniformly distributed in (0, 27). 
The signal received by the transceiver 


y(t) = g2e(2) 
has a power level in dBm of 





A? 2,2 
y = 10 log ae (1 + b? — cos ¢) 
=X, + X2+ X3+ UV, (4) 
where 
w = 10 log(1 + b? — 2b cos ¢) : (5) 


is the contribution of the echo canceler to the receive level of the VPA 
tone. Its effect can be evaluated by deriving the distribution of w. 
Since the density functions of 6 and ¢ are known, the density function 
of w can be shown to be 
fui(w), 10 log 0.25 < w < 10 log 0.75, 
fu(w) = fu2(w), 10 log 0.75 < w <0, (6) 
fus(w), 0 < w < 10 log 2.25, 


where : 


furi(w) = ky | as dv, 


1.25-q, 7492 





1 1.25—q3 
fuo(w) = hy | 1595 oy + | 96 yy 





Jimq, 4492 Tq, 4H 
1.25—-q3 
fus(w) = hy i he dv, 
~1 qq 
i! 
kh =——-,, 
V2nnog 


gqi=Ut+ vv? — 1+ qa, 
Q2 =u— vu? —-1+4q3, 


gs = 10", 
gs = V(v? — 1+ qs)(1 — v?), 
gs = exp[—(20 log g2 + 1)*/20"], (7) 
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and 

ge = exp[—(20 log qi + p)?/207]. 
The density function f,,(w) is plotted in Fig. 4. Since the mean value of 
1 + 6? — 2b cos ¢ is slightly larger than one, the expected value of w 
is slightly larger than zero. 

The receive level y in eq. (4) is the sum of random variables with 
known density functions; therefore, its distribution is also known. The 
transceiver detector level u is assumed to be uniformly distributed in 
(—9.4, —6.6) dBm. Thus, the probability of vpa failure with the echo 
canceler can be derived as 


Pue = P{y <u} 


10log 0.75 2—-w—-6.6 1 
-s ci [ ard | BI ty dw 
-« \J1010g0.25 w-9.4 1.25-q, 7442 
0 2—-w—6.6 939 1.25-q3 
3 
+ i ~ qQz ax (f = dv + i ae ao) dw 
1010g0.75 J z—w—9 Vi-q, 14% vixq, 9401 


10log2.25 rz—w—-6.6 1.25—g3 
ifs i | q7 dx | aids io aw} dz. (8) 
0 z-w—9.4 1 G1 


The probability of vpa failure without the echo canceler can be 
calculated as 








z—-6.6 
P, = P{x, + x2 + x3 <u} = —— qi dx dz, (9) 
TE ~o + 2—-9,4 
where 
xt+i1 x-l1 
we x ( 2p ) - e( 2p 
and 
V2x 
erf (x) = e?? dy 
a 


is the error function. Equations (8) and (9) are used to plot Fig. 5. 


APPENDIX B 
The probability of sTL is derived in this Appendix. An STL is 
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generated on a trunk if the trunk fails a vpa check, if another trunk in 
the same carrier group passes a subsequent vPA check, and if the first 
trunk fails a VPA recheck. For simplicity, the second trunk is assumed 
to be independent of the first trunk. The probability that a trunk fails 
two consecutive VPA checks will be derived first. It is assumed that the 
small probability of selecting the same transceiver for the two VPA 
checks is negligible. 

A vpa check on a trunk fails if x; + x2 + x3 + w <u. The values of 
Xe, X3, and w remains unchanged for a vPA recheck on the same trunk. 
Thus, the probability of interest is P{x2 + x3 + w <u — x}. Let z= 
Uu — X1, § = X2 + x3 and y = X2 + X3 + w. It can be shown that 








1 
56 (10.4 + 2), —8.4>2z2=-10.4, 
f-(z) : 76222 -8.4 
z = = 7? —(0=Z2= —-6.4, 
5.6 
J ( 5.6) 5.6 =>2z= —-7.6 
60000 Pe eae ee ee 
1 2/42 
f(s) = exp(—s*/4p’), (10) 
2p 
and 
10log0.75 fl 
Aopen, | | 29095 oy 
10log0.25 J125-q, 9442 
0 1 1.25—@3 
- b| {| Ot ay + | ae ao} 
1010g0.75 \J Vi—q, 94Q2 Nar q441 
10log2.25 1.25—q3 
=e | | 159895 a, ay 
0 ai G41 
where 
1 
| ee 
222" opg 
and 


gs = exp[—(y — w)?/4p’]. 


For a given value of y, a VPA failure is generated with the probability 
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Ply <2} -| f(z) dz 
: y 
1, y = -10.4, 
as («08 + 10.4y + *), —8.4>y => -10.4, 
5.6 2). 
1 
={ — 56 (13.2 + 2y), —76>y= —-8.4, 
es (15.06 + 5.6y + *), -5.6=>y=> —-7.6, 
5.6 2 
0, y= —5.6. 


Therefore, the probability of two consecutive vPaA failures on a trunk 
‘with echo canceler is 


Pr -| Ply <2}P{y<z}f(y) dy 


-10.4 1 (7 2\ 2 
= | [Ay) dy += [ («8.48 + 10.4y + *) 
_ . -10.4 


—-7.6 


1 
x f(y) dy + =a (13.2 + 2y)*f,(y) dy 
. —8.4 


1 —5.6 y? 2 


—71.6 
The probability of stL with echo canceler is then 
Pre = Pre (1 ba: Pe): (12) 


Without echo canceler the probability of two consecutive vPA failures 
on a trunk, Py», is obtained by replacing f,(y) in eq. (11) by f:( y) in (10). 
The probability of stL without echo canceler is then 


P, = P2-(1 — P.). (13) 
Equations (12) and (13) are used to plot Fig. 6. 


APPENDIX C 


The probability of scGa is studied in this appendix. An SCGA is 
generated if two trunks in the same carrier group fail successive VPA 
checks. If the two trunks are independent, the probability of scGa is 
simply the square of the probability of vpA failure. This is plotted as 


ECHO CANCELER 1329 


lower bounds in Fig. 7. Since the two trunks share identical carrier 
facilities from the switching offices to the earth stations, they may be 
dependent in some unknown way. A pessimistic estimate of the prob- 
ability of SCGA is to assume that x2 and x3 remain the same for the two 
trunks selected. Thus, the probability of itnerest is P{x; + w—u< 


— (x2 + x3)}. 
Let ¢= x, + wandh=t-— u. Then, 


fa(d), 10 log 2.25 ~ 1<¢< 10 log 2.25 + 1, 
fe(d), 1<¢< 10 log 2.25 — 1, 


fis(t), 10 log 0.75 +1<t<1, 


f(t) =} fra(t), —1<t< 10 log 0.75 + 1, 


fist), 10 log 0.75 -1<¢<-1, 


fis(t), 10 log 0.25 + 1<¢< 10 log 0.75 — 1, 
felt), 10 log 0.25 — 1 << 10 log 0.25 + 1, 


where 


1 10log 0.25 
fa(t) = | fus(w) dw, 
t-1 
t+1 


f(t) =5 Oe 


t-1 


1 t+1 1 0 
fs(t) = | fws(W) dw + | fw2(w) dw, 


1 é+1 
fault) = 5 [ foa(w) dw + ; | 


101log 0.75 
10log 0.75 


fs(t) -| fu2(w) dw + | fur(w) dw, 


0log0.75 t-1 


0 


f(t) = | fui(w) dw, 
t-1 
and 


f(t) = 5 i Palwidin: 


Olog 0.25 


The density function of h is given by 


10log 0.75 


fu2(w) dw + | fui(w) dw, 


t-1 
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frlh) = 


where 


fal), 10 log 2.25+ 84<h< 10 log 2.25 + 10.4, 
fn2(h), 10 log 2.25+ 7.6<h< 10 log 2.25 + 8.4, 
fas(h), 10.4 <h < 10 log 2.25 + 7.6, 
fra(h), 10 log 0.75 + 10.4<h < 10.4, 

frs(h), 10 log 0.75 + 5.6< h < 10 log 0.75 + 10.4, 
fas(h), 8.4 < h < 10 log 2.25 + 5.6, 
fiilh), 76<h< 8.44, 

fas(h), 10 log 0.75 + 84<h< 7.6, 

fao(h), 10 log 0.75 + 7.6<A< 10 log 0.75 + 8.4, 
frrolA), 5.6 < h < 10 log 0.75 + 7.6, 
faulh), 10 log 0.25 + 10.4<h< 56, 

fre(h), 10 log 0.75 + 5.6< A < 10 log 0.25 + 10.4, 
fnis(h), 10 log 0.25+ 84<h< 10 log 0.75 + 5.6, 
fnia(A), 10 log 0.25 + 7.6<h < 10 log 0.25 + 8.4, 
fas(A), 10 log 0.25+ 5.6<h< 10 log 0.25 + 7.6, 


10 log 2.25+1 
fulh) =s3 | Iu(t) dt, 
2.8 h-9.4 


l 10log 2.25+1 l 10log2.25—1 

fao(h) = | fu(t) dt +38 f(t) dt, 
“ J 10log 2.25—1  Sh-9.4 
1 h-6.6 1 10log 2.25—1 

fas(h) = | fa(t) dt | fa(t) dt, 
° 10log 2.25—1 ° h-9.4 


h-6.6 10 log 2.25—1 
fralh) = 55 | fu(t) dt + f(t) at 
2.8 10log 2.25—1 2.8 1 


1 1 
Bema fis(t) dt, 
2.8 x 


h-6.6 10 log 2.25—1 


1 
fas(h) = 28 fault) dé +38 fio(t) dt 


; 10 log 2.25—1 
1 1 10log0.75+1 
+ | fis(t) dt + — fea(t) dt, 


2.8 10log0.75+1 2.8 h-9.4 
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For a given value of s, a vPA failure is generated with the probability 
P{h <-—-s} = | fr(h) dh 


0, for —s < 10 log 0.25 + 5.6, 


[. fais(h) dh, 
10log 0.25+5.6 


for 10 log 0.25 + 5.6 < —s < 10 log 0.25 + 7.6, 


10 log 0.25+7.6 —s 
| fas(h) dh + { fris(h) dh, 
1 


Olog 0.25+5.6 10 log 0.25+7.6 


for 10 log 0.25 + 7.6 < —s < 10 log 0.25 + 8.4, 


101log0.25+7.6 10 log 0.25+8.4 
| fus(h) dh +| fris(h) dh 


10log 0.25+5.6 10log 0.25+7.6 


+ [ fnsth) dh, 
10log0.25+8.4 


for 10 log 0.25 + 8.4 < —s < 10 log 0.75 + 5.6, 


10log0.75+5.6 
[ tos | tau | fris(h) dh 
1 


Olog 0.25+8.4 


+ | fri2(h) dh, 
10log0.75+5.6 


for 10 log 0.75 + 5.6 < —s < 10 log 0.25 + 10.4 


10 log 2.25+8.4 
| fo + | frat: +f faa(h) dh 
10 log 2.25+7.6 


+ i fulh) dh, 
10 log 2.25+8.4 


for 10 log 2.25 + 8.4 < —s < 10 log 2.25 + 10.4 


1, for 10 log 2.25 + 10.4 < —s. 
Let p = —s. The pessimistic estimate of the probability of scGa with 
echo canceler is given by 
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| P{h < p}P{h <p} fs(p) dp 


1010g0.25+7.6 / fp 2 
7 { ( | fns(h) at) fap) dp 
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Olog 0.25+5.6 10log 0.25+5.6 


1010g0.25+8.4 101og0.25+7.6 
+ | ( | fais(h) dh 
1 1 


Olog0.25+7.6 Olog 0.25+5.6 


. 2 
a | fra) an) fp) dp 


10log 0.25+7.6 


10 log 2.25+10.4 
et] (| fs | fit + 
10 log 2.25+8.4 


10 log 2.25+8.4 p 2 
“3 [ fro(h) dh+ i filh) an) fs(p) dp 


Olog 2.25+7.6 10log 2.25+8,4 


+ i (5.6)"f.(p) dp. (14) 
1 


Olog 2.25+10.4 


Without echo canceler, the pessimistic estimate of the probability of 
SCGA is simply P;, the probability of sTL given in eq. (13). Equations 
(13) and (14) are used to plot the upper bounds of Fig. 7. 
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Profile Characterization of Optical 
Fibers—A Comparative Study 
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The refractive index profiles of several multimode optical fibers 
were measured by four of the current state-of-the-art techniques. 
These include interference microscopy on slab samples, transverse 
interference microscopy on whole fiber samples, the focusing method, 
and the refracted near-field method. The profile of the parent preform 
of one of the fibers was also measured by both the focusing method 
and by a ray-tracing approach. Comparisons of the results and the 
measurement methods are made emphasizing the applicability and 
use of the various techniques. 


1. INTRODUCTION 


Precise methods for measuring index profiles in both multimode and 
single-mode optical fibers and preforms are required if the desired 
ideal index distributions are to be produced. This paper will provide a 
comparison of four current state-of-the-art techniques for making 
these measurements on fibers and two similarly current methods used 
for preforms. Implications of the results especially related to the 
applicability and utilization of the various techniques will be discussed. 

It is important to have some feeling for the structure of the object 
being profiled. The fibers and the preform studied were produced by 
the modified chemical vapor deposition process.’ In this procedure a 
silica tube is mounted on a glass-working lathe and slowly rotated 
while reactants and dopants flow through it in an oxygen stream. An 
oxy-hydrogen burner is slowly traversed along the outside of the tube 
to provide simultaneous deposition and fusion of a layer of the reacting 
materials. On the order of 50 layers are deposited by multiple passes 
of the burner. To fabricate graded-index fibers, the dopant concentra- 
tion is gradually increased with increasing layer number. 

At the conclusion of deposition, the temperature of the burner is 
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raised to collapse the tube into a solid preform. A diagnostigram’ of 
the remaining end of the preform after being pulled into a fiber is 
shown in Fig. la. This display is obtained by expanding and collimating 
the light from a CW laser and allowing it to impinge upon the preform. 
The diagnostigram provides a simple and rapid nondestructive means 
of investigating internal layer structure. Details of its operation can be 
found in the Appendix. A shadow of the blunted end of the preform is 
seen at the lower right and each of the deposition layers is seen as a 
horizontal line. The uppermost bright line is the core-cladding bound- 
ary and the region immediately below it is a 320-ym-thick barrier layer. 
The core radius is about 3.5 mm and the length of preform seen in the 
diagnostigram is 12 cm. The point to note is the rich structure and 
resolvability of the deposition layers. Another view of this structure 
can be obtained by immersing the preform in index-matching oil and 
observing the incoherently illuminated core with a video camera, as is 
done in the focusing method.’ This representation, shown in Fig. 1b 
for a several millimeter length of preform, also emphasizes the individ- 
ual layer structure, which is in addition displayed by the curve at the 
side of the display. 

These structural variations also appear, appropriately scaled, in the 
fiber in a one-to-one correspondence‘ and corroborates the fact that 
the same distribution of refractive index that is introduced into the 
preform exists in the fiber.’ Generally, the scale of the variations in 
the fiber is on the order of less than a wavelength and they are, 
therefore, not observed by the measurement technique. Notable ex- 
ceptions occur near the axis where the deposition layers are thickest 
and in any region where either several layers have the same index or 
thicker than normal layers are produced, due to fabrication faults. 

The profile, built up by varying the dopant concentration in each of 
the layers, can be measured by a variety of techniques. How accurately 
the measurements should be made and which technique is best for 
making them are always difficult questions to answer since they 
involve trade-offs of many factors. To be included are time, cost, and 
effort considerations, all of which one would like to minimize, and 
sensitivity, accuracy, and resolution which one would like to maximize. 

A handle on the question of accuracy is provided by the following 
considerations. The theoretical bandwidth that can be realized with 
an optimum profile is about 8000 MHz x km for a fiber with a 
maximum index difference value of 0.02.° To achieve this high band- 
width requires that the exponent, g, of the power-law profile have a 
definite optimum value near g = 2 (for germanium dopant and an 
operating wavelength of 0.9 wm). A departure of only 0.05 from this 
optimum g value is sufficient to reduce the fiber bandwidth by more 
than one order of magnitude. Clearly the profiling technique must 
determine g to better than 0.05 if a meaningful correlation between 
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Fig. 1—Layer structure of mcvpD fabricated preform as observed by: (a) diagnostri- 
gram, and (b) index immersion. 
fiber performance and index profile is to be obtained. In order to 
achieve this accuracy, An(r) [the difference between the refractive 
index of the fiber core and its cladding value] must be measured with 
a precision of 1 part in 10%. 

Very slight local distortions of the refractive index profile from its 
optimum shape also decrease the fiber bandwidth markedly. A distor- 
tion of ten sinusoidal periods over the fiber radius with an rms 
amplitude of 0.6 percent of the maximum index difference reduces the 
bandwidth from 8000 MHz x km to about 200 MHz x km.®° An rms 
distortion amplitude of only 0.15 percent for the same ripples reduces 
the bandwidth to about 800 MHz x km. The precision of the An(r) 
measurement would again have to be about 1 part in 10* to detect even 
the 0.6 percent distortion. In addition, of course, the spatial resolution 
of the technique must be sufficient to resolve the perturbations. 


OPTICAL FIBER PROFILES 1337 


SLAB 






_ OUTPUT 
LIGHT 


Bee ae 
~MIRRORS ~ 


TRANS 






__ OUTPUT 
LIGHT 


NY 7 
~MIRRORS~~ 
















/MICROSCOPE , , MICROSCOPE, 
/ \ 7 N 


7 REFERENCE 


REFERENCE rae 
SLAB / 


~ ~ \ SEMITRANSPARENT 
MIRRORS 


~ — Ss SEMITRANSPARENT 
MIRRORS 


__INCIDENT 
LIGHT 


__ INCIDENT 
LIGHT 


(a) (b) 


Fig. 2—Schematic diagram of (a) slab interferometric method, and (b) transverse 
interferometric method. The interference microscope used in both cases is identical. 


ll. PROFILE MEASUREMENT METHODS 


The specific profiling methods used in this study are shown diagram- 
matically in Figs. 2 and 3. They are interference microscopy on SLAB 
samples (Fig. 2a); transverse interferometry on whole fiber samples 
(Fig. 2b); the focusing method (Fig. 3a); and the refracted near-field 
method (Fig. 3b). They are abbreviated by the terms SLAB, TRANS, 
FOCUS, and RNF, respectively. These methods will be discussed briefly; 
further details on their practical implementation can be found in the 
Appendix. 


2.1 Slab interferometry 


Interference microscopy on SLAB samples, utilizing the potential 
accuracy of interferometry was historically the first of these methods 
to be used,’ and is generally accepted as the method to which newly 
developed techniques are compared. The SLAB sample is cut from an 
encapsulated fiber (or preform tip) and polished so that the faces are 
flat and parallel. The cutting and polishing procedures are both diffi- 
cult and time-consuming. Special techniques are required to avoid 
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Fig. 3—Principle of operation of (a) focusing method, and (b) refracted near-field 
(RNF) technique. 





composition-dependent thickness variations, which can lead to sub- 
stantial errors.’ It is also necessary for the sample to be thin enough so 
that rays traversing it are not bent and focused, producing curved 
wavefronts and, hence, erroneous results. Sample preparation requires 
about one day. This time can be reduced on a per-sample basis by 
processing several different fibers, which have been epoxied into one 
capillary tube, at the same time. It should also be noted that this 
procedure is not inherently nondestructive. 

To observe the samples, of course, an interference microscope is 
required. Interference lens attachments to ordinary microscopes gen- 
erally involve passing the light through the specimen twice, thus 
compounding possible errors. Best results are obtained with a single- 
pass Mach-Zehnder geometry, but microscope cost and availability are 
additional major considerations in adopting this method. 

Relative index values accurate to about 2 parts in 10‘ can be realized 
routinely and by electronically processing the output of the microscope 
measurements relatively accurate to about 1 part in 10°, as necessary, 
for example, in profile dispersion work, have been achieved.’ Automatic 
computer processing of the output also serves to reduce analysis time. 
Spatial index resolution is somewhat limited in that it is not possible 
to combine maximum lateral resolution and exact phase measurements 
in a single instrument.’ 


2.2 Transverse interferometry 


Sample preparation can be eliminated by using the transverse inter- 
ferometric method (Fig. 2b). In this technique, a length of fiber is 
immersed in matching oil on the stage of the interference microscope 
and illuminated at right angles transverse to its axis. The matching oil 
removes the influence of the outer cladding boundary. The total optical 
path length of a light ray is expressed as an integral and the index 
distribution is obtained from the measured fringe shift by solving an 
integral equation.’ Unlike the sLaB approach, in which the entire core 
is accessible, transverse interferometry assumes circular symmetry 
and, hence, geometry variations, which can adversely affect the profile, 
will not be detected unless special care is taken to make several 
measurements for different rotational positions of the fiber. By auto- 
mating the measurement, index profiles of a fiber can be obtained 
within a few minutes after its manufacture. The accuracy of the 
method is about an order of magnitude less than that of the SLAB 
approach, and it is subject to a large error in the region near the fiber 
axis. On the other hand, this technique resolves detail in the fiber 
structure with higher resolution than the SLAB method, as will be seen 
later in the actual profiles. 
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2.3 Focusing method 


The focusing method” (Fig. 3a) is similar to transverse interfero- 
metry in that it is also nondestructive and uses transverse illumination. 
Otherwise they are very different. The focusing method does not 
require an interference microscope or rely in any way on interferome- 
try. Moreover, the technique is readily applicable with high accuracy 
to fiber preforms.* 

In this method, the fiber, observed with a microscope and the 
preform with a camera lens, are immersed in index matching fluid. The 
core, acting as a cylindrical lens, focuses the light, whose power density 
distribution in the observation plane is detected by a video camera. 
After digitizing the power density a computer calculates the index 
profile by solving an integral equation. The profiles so obtained from 
circularly symmetric cores are comparable in accuracy and resolution 
to those produced by interference microscopy of SLAB samples.” 

Extreme experimental care, however, is involved in the focusing 
method since it measures absolute light intensities. The optics, match- 
ing oil, and fiber (or preform) itself must be very clean; the video 
detector must be linearized and the incident light intensity must be 
uniform, either intrinsically or through a calibration procedure. 


2.4 Ray tracing 


A related method to measure the profile in preforms is ray tracing.“ 
Further details of this technique are also given in the Appendix. This 
method involves scanning a focused laser beam perpendicular to the 
axis of the index-matched preform and recording the exit angle of the 
beam as a function of distance from the axis. The profile of the preform 
is then reconstructed by taking the inverse Abel transform of the 
deflection function. Indeed, the mathematics of this and the focusing 
method are nearly identical’® and lead to very similar profiles, if an 
equal number of data points are measured and processed. 


2.5 Refracted near-field 


The refracted near-field (RNF) method” (Fig. 3b) relies on the power 
escaping sideways from the core into the cladding to determine the 
refractive index profile of the fiber. The fiber, immersed in a matching 
oil whose index of refraction is greater than that of the cladding, is 
passed through a small hole in an opaque disc. Part of the light focused 
into the fiber is guided while the rest appears outside of the fiber as a 
hollow cone. If all of the leaky modes, contained in the inner part of 
this cone are blocked by the disc, then the light passed varies linearly 
with the index of refraction of the fiber at the point at which the 
incident light is focused. Thus, by scanning the incident light across 
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Fig. 4—Two fiber-SLAB samples observed by (a,c) ordinary microscopy, and by (b,d) 
interference microscopy. 


the end of the fiber, the profile is obtained directly from the output of 
a detector that collects the light passing the disc.’° This method, as 
will be seen in the profiles, has excellent spatial resolution and also 
possesses the ability to analyze noncircular cores. The precision of the 
index measurement is estimated in one recent embodiment”® to be 
about 4 x 10°. 


lil, MEASUREMENT RESULTS 


The fiber samples used in this study were specifically chosen to 
possess a variety of features, which would severely test the limits of 
the different profiling methods. 

A SLAB sample of the first fiber is shown in Figs. 4a and 4b as 
observed with ordinary and interference microscopy, respectively. The 
fiber is seen to possess severe perturbation in layer structure through- 
out the core with marked variations especially near the center. The 
index profiles as obtained by the SLAB, FOCUS, and TRANS methods are 
shown in Fig. 5 by the solid, dotted, and dashed curves, respectively. 
The index and core radius values, which are in very good agreement, 
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are those obtained by each of the measurements and no scaling was 
employed. Generally the profile shapes are similar, with the TRANS 
profile exhibiting more detailed structure of the perturbations. Slight 
differences are accounted for by the assumption of circular symmetry 
in the FOCUS and TRANS cases and the lack of such an assumption for 
the SLAB. 

A comparison of the RNF profile (solid curve) and the SLAB profile 
(broken curve) for this fiber is shown in Fig. 6. The SLAB, FOCUS, and 
TRANS profiles were all measured by the author at Bell Laboratories, 
Crawford Hill. The RNnF profiles were provided by Jeff Saunders at 
Bell Laboratories, Atlanta.’ The resolution of structure in the RNF 
profile is striking in comparison to the SLAB, which appears as a near 
average through the RNF results. The scale here and in subsequent 
profiles has been eliminated for clarity and ease of comparison. The 
RNF profile also shows a steeply rising region, and an index depression 
at the core-cladding interface, features not well resolved by the SLAB. 
The SLAB measurement is also not able to resolve the central depres- 
sion due to its steep gradient. 

A comparison of the RNF profile (solid curve) and the TRANS profile 
(broken curve) is shown in Fig. 7. The superior resolution of the TRANS 
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Fig. 5—Profiles of fiber shown in Figs. 4a and 4b (fiber no. 1) as obtained by SLAB 
(solid), Focus (dotted), and TRANS (broken) techniques. 
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Fig. 6—Comparison of SLAB (broken) and RNF (solid) profiles for fiber no. 1. 
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Fig. 7—Comparison of RNF (solid) and TRANS (broken) profiles for fiber no. 1. 


measurement to that of the previously shown SLAB is clearly seen in 
that now the curves are practically identical, except for the region near 


the center where RNF shows greater detail. 
The comparison of RNF (solid curve) and Focus (broken curve) 
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Se 


Fig. 8—Comparison of RNF (solid) and Focus (broken) profiles for fiber no. 1. The 
dotted curve is a portion of a focused profile obtained by focusing within the core. 





Fig. 9—Comparison of RNF (solid) and sLAB (broken) profiles for fiber shown in Figs. 
4c and 4d (fiber no. 2). 
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Fig. 10—Fiber no. 3 as observed by (a) ordinary microscopy, (b) slab interference 
microscopy, and (c) the focusing method. 





Fig. 11—Comparison of RNF (solid) and SLAB (broken) profiles for fiber no. 3. 


shown in Fig. 8 is similar to the RNF versus SLAB results. A second 
FOCUS profile, a portion of which is shown by the dotted curve, was 
obtained by focusing within the core to obtain greater resolution. 
Indeed the resolution in this region is now comparable to RNF but the 
remainder of the profile (not shown) is distorted because of the 
violation of the focusing condition. Focusing within the core actually 
satisfies the focusing condition locally for those rapid variations that 
tend to focus the incident rays much more steeply than the remainder 
of the core. 

A second SLAB sample with somewhat smaller index variations is 
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shown in Figs. 4c and 4d, again by ordinary and interference micros- 
copy, respectively. A comparison of the RNF (solid curve) and SLAB 
(broken curve) profiles is presented in Fig. 9. Because of the smaller 
index fluctuations, the curves agree quite well but, again, the high 
resolution of RNF is apparent. 

A third SLAB sample of a fiber fabricated with several (~10) discrete 
changes in dopant concentration as a function of core radius is shown 
in Figs. 10a and 10b by normal and interference microscopy, respec- 
tively. Figure 10c shows the whole fiber sample as observed by the 
focusing method. 

A comparison of the RNF (solid curve) and SLAB (broken curve) 
profiles seen in Fig. 11 shows good agreement as far as the general 
profile shape is concerned, but the SLAB curve is, again, almost an 
average through the well-resolved layer structure of RNF. The FOCUS 
profile is similar to the SLAB, except when it is obtained by focusing 
within the core in which case specific local features can be brought 
out, at the expense of an overall distortion, with resolution comparable 
to RNF. 

The TRANS profile of this fiber, shown by the broken curve in Fig. 
12, is in excellent agreement with the RNF profile (solid curve). The 
resolution of the layers, the widths of which are somewhat larger than 





Fig. 12—Comparison of RNF (solid) and TRANS (broken) profiles for fiber no. 3. Note 
the comparable resolution. 
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Fig. 13—Fiber no. 4, possessing a strong index perturbation, as observed by the 
various measurement methods. 


in the first sample, is equally good. The index depression and steep 
rise in the index distribution are also the same in both. 

A fourth sample having a major index perturbation at 50 percent of 
the radius is shown in Figs. 13a and 13d by ordinary microscopy, SLAB 
interferometry, focusing method, and transverse interferometry, re- 
spectively. The perturbation is clearly seen in each method of obser- 
vation. A comparison of the SLAB (broken curve) and RNF profile (solid 
curve) is shown in Fig. 14. The profile shapes are similar but the index 
distortion is more prominent in the RNF measurement. The TRANS 
profile, shown as the broken curve in Fig. 15, on the other hand, shows 
the perturbation as clearly resolved as RNF (solid curve). The TRANS 
profile also displays the fine index variations, as does the RNF result, 
which lie closer to the center of the core. These fluctuations are absent 
in the SLAB and also in the Focus profiles both of which are very 
similar. 

Finally, a fiber with a relatively perturbation-free index distribution 
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Fig. 15—Comparison of RNF (solid) and TRANS (broken) profiles for fiber no. 4. 


was studied. The fiber as viewed by the SLAB, TRANS, and FOCUS 
techniques is shown in Fig. 16. To be noted in particular is the very 
uniform appearance (except for dirt specks) of all the samples indicat- 
ing a lack of index distortions. The profiles obtained by the different 
methods, as might be imagined, are all very similar. The RNF (solid 
curve) and SLAB profiles (broken curve) are shown in Fig. 17 and the 
RNF (solid), Focus (dotted), and TRANS (broken) profiles are seen in 
Fig. 18. The only differences in the index distributions is a steep initial 
rise in the profile and some very slight ripple seen in the RNF result 
near the core center. This particular fiber also had a strong asymmetry 
in the index distribution near the center. This is not apparent in the 
FOCUS and TRANS results since they depend upon the assumption of 
circular symmetry. 
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Fig. 16—Fiber no. 5 possessing a relatively perturbation-free profile as observed by 
the various measurement methods. 





Fig. 17—Comparison of RNF (solid) and SLAB (broken) profile for fiber no. 5. 
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Fig. 18—Profiles of fiber no. 5 as measured by RNF (solid), Focus (dotted), and TRANS 
(broken) techniques. 


An important question is the relationship of profiles measured in 
the fibers to the actual profile existing in the corresponding preforms. 
A knowledge of this relationship would shed additional light on the 
applicability of the various methods and give confidence that structural 
features observed in the fiber profiles are not measurement artifacts. 
Of particular interest is the steep initial rise in the index distribution 
of the last fiber as seen by RNF but not clearly defined with the other 
methods. Is this feature in the preform or not? 

To answer this question the preform corresponding to the last 
sample was profiled by the two methods previously described, the 
focusing method as applied to preforms’ and the ray-tracing technique 
as used at Western Electric’s Engineering Research Center."* 

The preform profiles are shown in Fig. 19; the solid curve is obtained 
by the ray-tracing technique and the broken curve by the focusing 
method. The profiles can barely be distinguished, except in the region 
near the center where the previously mentioned index perturbation 
makes the profile shape orientation dependent, and in the resolution 
of the finer deposition layers, which the ray-tracing method achieves 
by processing about ten times as many data points as the focusing 
method. 

A comparison of the focused preform profile and the RNF fiber 
profile is shown in Fig. 20. Scaling of the radial coordinate, of course, 
was performed but not of the index values. The agreement of the 
profiles is excellent. It is seen that the initial steep rise which appears 
in the RNF profile is indeed present in the preform. It does, however, 
only appear as a single step in the RNF result, whereas in the preform 
it has a fine double-step structure. The chosen resolution of the focused 
preform profile is, thus, seen to correspond to the actual fiber profile. 
The greater detail included in the ray-tracing technique is absent even 
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Fig. 19—Profiles of the preform stub from which fiber no. 5 was pulled as measured 
by the ray-tracing technique (solid curve) and the focusing method (broken). The two 
halves of each profile were obtained from measurements at two orientations. 


from the RNF-profiled fiber, because of its small-scale structure. In 
general, a very high degree of correspondence between fiber and 
preform profiles exists. 

The various fiber profiles were further compared by taking the rms 
deviation of their differences as a function of radial position. This is 
the severest possible comparison in that absolute point-by-point values 
are looked at. Thus, a slight shift of a feature could result in a sizable 
difference. A second less restrictive comparison was also made by 
computing the best-fit g values to the various profiles. This emphasizes 
the profile shape at the expense of the location and existence of small 
perturbations. 

It was found that the various profiling methods give just about 
identical results for relatively smooth profiles. The rms deviations of 
the profiles of the last fiber (shown in Figs. 16-20), were all less than 
one percent, excluding the index depression in the center. The central 
depression gives rise to a few percent difference on its own because of 
the different ways it is resolved by the various methods. The best fit 
g curves were all within 0.05 of each other. The RNF profile has a g 
value of 2.024 and that of the preform (as measured by the focusing 
method) a value of 1.990, a difference of less than 0.035. 

As expected, the rms differences between the profiles increase for 
those possessing rapid variations. These differences for the first four 
fiber samples shown amounted to three to five percent, again excluding 
the index depression. The g values of the respective profiles were all 
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Fig. 20—Comparison of preform profile (broken) and corresponding RNF fiber profile 
(solid) for fiber no. 5. 


within a few percent, as would be expected from the similarities of the 
curves. 

We conclude then, that as long as the profile is smooth all of the 
measurement methods are equivalent. For profiles containing rapid 
perturbations, the RNF and TRANS methods are about comparable with 
the former doing better in resolving very rapid variations. The RNF 
method also has the advantage of not requiring an interference micro- 
scope nor elaborate computer processing of the data. On the other 
hand, it does require a separate calibration procedure and reasonable 
care with cleanliness of the optics. 

The fact that the rms deviations are less than one percent and the 
best-fit g values are within 0.05 for relatively perturbation free profiles 
lends confidence to the ability to make valid bandwidth predictions 
for current fibers of this type based on the various measured profiles."” 
Deviations of this magnitude, while reducing the bandwidth by about 
one order of magnitude, still result in bandwidths close to 1 GHz. As 
profiles get even smoother these deviations will presumably decrease 
and allow for even better predictions. Meanwhile, improvements in the 
accuracy of the profiling methods themselves, can also serve to reduce 
these interprofile deviations and lead to even better agreement. An 
improvement in accuracy of about a factor of 5 is required to mean- 
ingfully measure deviations of the ideal profile. On the other hand, 
one then enters the realm of theoretical uncertainty as to what the 
ideal ideally is. The current state-of-the-art of index measurements 
should then be able to go a long way in providing the feedback 
necessary to improve current profiles before their limitations are 
indeed felt. 

While we have presented four current profiling techniques, there 
are, of course, other existing methods each subject in use to their own 
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set of trade-offs. New methods both for fibers and preforms are 
reported regularly and undoubtedly as they prove their value will find 
use in the important task of index profiling. 
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APPENDIX 
Implementation of Profile Measurement Methods 
Preform diagnostigrams 


Preform diagnostigrams provide a sensitive, nondestructive, and 
noncontacting means of obtaining structural information. This infor- 
mation includes a measurement of core size and core eccentricity, core 
cladding interface structure, individual deposition layer structure and 
variations, imperfections within the core and the cladding, and the 
presence of an axial refractive-index depression. 
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Fig. 21—Arrangement to produce preform diagnostigrams. 
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A diagnostigram is produced by the arrangement shown in Fig. 21. 
Light from a He-Ne laser is incident upon an oscillating mirror that 
serves to transform it into a line. The line of light is expanded and 
collimated by two cylindrical lenses and traverses the preform. The 
beam is about equal in width to the diameter of the preform, and its 
length can be varied as desired by adjusting the amplitude of the 
oscillating mirror. Typically a length of 15 cm is used. An arrangement 
of cylindrical lenses replacing the oscillating mirror can also be utilized. 
The light traversing the preform is then incident upon either an 
observation screen or photographic film. 

The pattern, as shown in Fig. la (of the text), consists of bright and 
dark lines. The width of the bright lines represents the actual geometric 
width of discontinuous index steps in the core, while the width of the 
dark lines represents the amount An of refractive-index discontinuity. 
The discontinuities arise during the deposition process and represent 
distinct deposition layers. There exists a one-to-one correspondence 
between the observed lines and the deposited layers. Further details 
can be found in Ref. 2. 


Slab and transverse interferometry 


In slab interferometry the fiber sample is placed in one arm of the 
interference microscope and a homogeneous reference SLAB with re- 
fractive index nz is placed in the reference arm (Fig. 2a). An arrange- 
ment of practical implementation is shown in Fig. 22. Figure 23 displays 
the fringe shifts observed in a graded-index sample, the shift S of a 
fringe depending on its position in the fiber core, S = S(r). The 
difference between the refractive indices of core and cladding can be 
expressed in terms of this fringe shift S(7), the uniform fringe spacing 
in the cladding D, the vacuum wavelength of light A, and the SLAB 
thickness ¢ as 

AS(r) 


n(r) — ne = Dr (1) 





To measure the fringe shift a video camera looks into the interfer- 
ence microscope and sends its electrical output signal to a video 
digitizer. The 8-bit digitizer is computer controlled and addresses 
specific, preselected points in the video field. A video monitor and a 
plotter for recording the processed information—the index profile, is 
also included. The computer directs the vertical sample line seen in 
Fig. 23 to collect information on either side of the core (along lines A- 
B, and C-D), which is then used to determine the fringe spacing and 
to compensate for a tilt of the entire fringe pattern. The computer then 
advances the sample line in small increments, moving it through the 
core region, measuring the displacement of the fringe that goes through 
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Fig. 22—Arrangement for automatic refractive index profiling of SLAB samples. 





Fig. 23—Fringes observed in a graded-index SLAB sample. 
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Fig. 24—Fringes observed in a graded index whole fiber sample observed transversely. 


the core center. The computer determines the fringe positions by 
counting and searching for the minimum light level whose location it 
pinpoints by least mean square fitting of a parabola using a number of 
points in the vicinity of the minimum. 

The fringe displacement is recorded as a function of the radial 
coordinate 7 measured from the core center and the resulting function 
S(r) is used to compute n(r) — nz according to eq. (1). The index 
distribution is then sent to the plotter. 

In transverse interferometry the fiber, after being stripped of any 
jacket, is again placed in one arm of the interference microscope (Fig. 
2b in text). The fiber is covered with a drop of matching oil into which 
the microscope objective is dipped. The reference branch contains only 
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a drop of matching oil. Each light ray incident upon the sample now 
passes through regions of varying refractive index and the total path 
length must be expressed as an integral. The refractive index difference 
between core and cladding is given by 








_. _ A [* |dS(e)| dp 
n(r) m= 3 | dae veor (2) 


in which p is the radial coordinate and rotational symmetry is assumed. 

The output field of the interference microscope now appears as in 
Fig. 24, and the fringe shift is measured with the computer-controlled 
video digitizer system just described. Processing of the fringe shift 
information requires that a numerical differentiation is performed first, 
followed by the numerical integration indicated by eq. (2). 


The focusing method 


The focusing method is shown in application to fibers and preforms 
in Fig. 25. The technique uses incoherent filtered light in transverse 
illumination. The fiber, observed with a microscope, and the preform, 
observed with a camera lens, are immersed in index matching fluid to 
avoid the deleterious influence of the outer cladding boundary. The 
core, acting as a cylindrical lens, focuses the light whose power density 
distribution in the observation plane is detected by a video camera. 
The observation plane is defined by the object plane on which the 
camera is focused. This plane must not be inside the fiber core, and it 
must not be placed so far away that rays have already crossed over 
after leaving the core. Good results are obtained when the observation 
plane is placed just outside of the core. 

The image of a preform seen by the camera is shown in Fig. 1b of 
the text. A monitor display of a fiber is shown in Fig. 26. 

The refractive index distribution is obtained by solving the integral . 
equation 
ma { t—yt) 


nr) — ne = 


| (3) 


The various parameters are defined in Fig. 3a of the text. 

The function y(t) is obtained from a measurement of the light power 
density distribution in the observation plane. This distribution is 
collected along the sample line and digitized under computer control 
as described previously. The computer also solves the integral equation 
and plots the resultant index profile. 


Refracted near-field method 


In this technique” (shown in Fig. 3b of the text) a light beam is 
focused on a spot at a distance r from the fiber axis with a convergence 
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Fig. 25—Arrangement for automatic refractive index profiling of fibers and preforms 
by the focusing method. 


angle that is larger than the acceptance angle of the fiber core. The 
light escaping from the core is partly contributed by power leakage 
from leaky modes. This part of the radiated power is blocked by a 
circular aperture which prevents light leaving below a minimum angle 
@min from reaching the detector. The refractive index difference be- 
tween the core and cladding, which is obtained from the light passed 
by the aperture, is given by 
n(r) — nz = n2Cos Ohrin[ COSmin — COS BJ (4) 
In this expression 9” refers to the input angle, P(r) is the light power 
reaching the detector, as a function of the position of the input beam 
and P(a) is obtained from the P(r) curve as the light power detected 
when the input beam is focused into the cladding. 
The experimental apparatus used in the implementation of this 
method by J. Saunders is shown in Fig. 27."° Light from a 5-mW He- 
Ne laser passes through a quarter wave plate and is focused onto a 50- 
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Fig. 26—Graded-index fiber as observed by the focusing method. 


pm pinhole. The light from the pinhole is then focused by a 20x 
microscope objective onto the end face of the fiber which is held in a 
moveable cell containing index immersion liquid. 

The disc that blocks the leaky rays is 1.38 cm in diameter and is 
supported by means of three fibers. The light passing the disc is 
directed by lenses to a large area detector whose output provides the 
profile. The microscope and TV camera provide a magnified view of 
the fiber core for alignment and monitoring purposes. 


Ray-tracing method 


The practical implementation of this method by L. S. Watkins, is 
shown in Fig. 28.’ A narrow beam from a He-Ne laser is reflected off 
a rotatable mirror and is focused by a lens through the index-matched 
preform. Rotating the mirror moves the ~20 pm beam across the 
preform. 

The deflected beam is collected by a lens and focused onto a linear 
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position sensor at its back focal plane. The output of the sensor is 
analyzed to give a voltage proportional to the deflection angle, which 
is then computer-processed in a similar manner to the focusing method 
to give the profile. 
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Fig. 27—Experimental apparatus for implementation of the RNF method. 
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Fig. 28—Arrangement for profiling of preforms by the ray tracing method. 
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Statistical Behavior of Crosstalk Power Sum 
With Dominant Components 


By S. H. LIN 
(Manuscript received February 4, 1981) 


The literature of digital transmission on wire-pair cables generally 
considers the probability distributions of both pair-to-pair crosstalk 
loss and its power sum to be normal on a dB scale. This paper 
presents extensive measured data of crosstalk among connectors and 
wire pairs on the backplane and associated stub cable of 466-type 
apparatus cases of the existing T1 system. The measured probability 
distribution of crosstalk power sum “bends” toward more severe 
crosstalk levels in the lower tail region, which is important for T1 
system engineering. This bend is because of the effects of a few 
dominant components (i.e., within-slot or within-harness crosstalk) 
in the power sum. The simple normal model is too optimistic by 4 dB 
in estimating apparatus-case-crosstalk power sum at 0.1 percentile 
level. This paper shows that both the Monte Carlo and the lower 
bound methods for power sum calculations predict this bend in close 
agreement with the measured data. Although apparatus-case- 
crosstalk power sum is worse than previously assumed, the perform- 
ance of T1 system has been adequately protected by the extra margin 
in the previous engineering rules to cover unknowns. 


l. INTRODUCTION 


Crosstalk interference is a prime limitation on the transmission 
capacity and the performance of digital transmission systems, such as 
T1,' T1C,”* SLC-40,* T1D,° and SLC-96,° on twisted multipair cables. 
An important step in the design of digital systems and their associated 
engineering rules is the characterization of the power sum of pair-to- 
pair crosstalk loss. The crosstalk power sum is the total crosstalk 
interference which appears on a given pair as a result of coupling from 
all disturbers on other pairs. 

The crosstalk power sum of a T-carrier system can be decomposed 
into two components: one component originates from crosstalk among 
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wire pairs in the cable, and the other component originates from 
crosstalk among wire pairs on the backplane of the repeater apparatus 
case and the associated connectors and stub cable. At each repeater 
location of a T1 system, an apparatus case is used to house 50 
regenerators of 50 one-way T1 systems. Figure 1 shows the 25-slot 
arrangement of the 466-type (without lightning protection device) 
apparatus case. Each slot holds two T1 regenerators. Figure 2 shows 
a portion of the wiring arrangement between the stub cable and the 
repeater connectors on the backplane of a 466-type apparatus case. 
The crosstalk originating from the wire pairs on the apparatus case 
backplane, connectors, and stub cable is known as apparatus-case- 
crosstalk (ACXT). 

In new, 800-series, plastic apparatus cases, a carefully controlled 
wiring layout is used to minimize the AcxT to such an extent that ACxT 
can be neglected in the T1 system engineering rules. However, exten- 
sive laboratory and field measurements indicate that the old vintage 
apparatus cases, such as 466-type, are major contributors to T1 inter- 
system crosstalk. In this paper, we study the AcxT data of 466-type 
apparatus cases in detail because this is one type of apparatus case 
which has been widely deployed in the existing T1 plant. The statistics 
of ACXT are, therefore, important in characterizing the performance of 
the existing T1 systems and future higher bit rate systems proposed to 
be used in the existing T1 environment. All the AcxT data presented 
in this paper were measured at 0.772 MHz, the Nyquist frequency of 
the T1 bit rate. As explained in Section II of Ref. 7, the extreme tail 
region (i.e., 0.1 to 0.025 percent) of the probability distribution of the 
crosstalk power sum is important in the engineering of digital trans- 
mission systems in twisted pair cables. 

Figure 3 shows the distribution of the power sum of AcxT of 466- 
type apparatus cases obtained from extensive laboratory® and field 
measurements.” In the simple normal model, the power sum data on 
Fig. 3 would be fitted by a straight line and the predicted 0.1 percentile 
would be 59 dB. However, the measured power sum distribution on 
Fig. 3 has a noticeable “bend” towards more severe crosstalk levels in 
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Fig. 1—The 466-type apparatus case with 25 repeater slots (ie., retainers) for T1 
repeaters. 
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Fig. 2—Wiring diagram on the backplane of 466-type Tl apparatus case showing 
within-slot crosstalk and non-within-slot crosstalk for slot 24. 


the lower tail region (<5 percent). It cannot be described precisely by 
any simple model such as normal or gamma. This paper shows that 
the bend is because of the dominant effect of the within-slot pair-to- 
pair crosstalk which is, on the average, 15 dB worse than the non- 
within-slot pair-to-pair crosstalk. It is demonstrated that both the 
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Fig. 3—Distribution of power sum of 50 pair-to-pair crosstalk losses of 466-type T1 
apparatus case from laboratory measurements in Atlanta and field measurements in 
Illinois, California, and Texas. 
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Monte Carlo and the lower bound method for power sum calculations 
predict this bend. Both methods and the data indicate that the 0.1 
percentile of the power sum distribution for 466-type cases is 55 dB 
which is 4 dB more pessimistic than the 59 dB predicted by the simple 
normal model. Therefore, the simple normal model may be too opti- 
mistic in estimating the Tl system margin. The previous engineering 
rules’” for T1 system contain extra margin to cover “unknowns.” The 
effect of ACxT on the bit-error-rate performance of T1 system has been 
adequately protected by the extra margin. The development of a more 
accurate ACXT model will reduce the unknowns and enable a greater 
exploitation of the system’s capability by mitigating the need for large 
“uncertainty” margins. M. H. Meyers" has also investigated acxT by 
a different approach. 


ll. ATLANTA DATA AND THE SIMPLE NORMAL AND GAMMA MODELS 


The AcxtT data of eight 466-type apparatus cases were measured in 
Bell Laboratories in Atlanta by using a computer operated transmis- 
sion measurement set.*)”'* 

The laboratory data are shown as circles on Figs. 4 and 5 for pair-to- 
pair crosstalk loss and the power sum, respectively. The solid line and 
the dashed line on Fig. 4 represent the gamma and normal approxi- 
mation, respectively, to the pair-to-pair distribution. The power sum 
distribution is strongly controlled by the behavior of pair-to-pair 
distribution in the tail region of low crosstalk loss.’ The gamma 
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Fig. 4—Distribution of pair-to-pair crosstalk loss of 466-type T1 apparatus case. Data 
measured from eight apparatus cases in Bell Laboratories, Atlanta. 
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Fig. 5—Distribution of power sum of 50 pair-to-pair crosstalk losses of 466-type T1 
apparatus case. Data measured from eight apparatus cases in Bell Laboratories, Atlanta. 


approximation fits the pair-to-pair data very closely in the critical 
region of low crosstalk loss, whereas the normal approximation is too 
pessimistic in this important tail region. 

The worst value of 60 dB on Fig. 4 does not imply that the pair-to- 
pair distribution is truncated at the 60-dB level. A finite sample 
measurement of a random variate (e.g., crosstalk loss) always yields a 
finite worst value even if the parent distribution of the variate is 
untruncated. The sample worst value varies randomly from one set of 
measurement (e.g., from one cable) to another. The probability distri- 
bution of the worst value (i.e., the extreme value) is the subject of 
extreme value statistics which have been studied extensively.'*’°"* 
Therefore, the existence of a finite worst value from a finite sample 
measurement of cable crosstalk loss does not necessarily imply that 
the parent distribution of crosstalk loss is truncated. 

Figure 5 shows that the power sum distribution predicted by the 
gamma model (solid line) agrees reasonably well with the measured 
data over a large portion of the distribution, but the discrepancy in the 
lower tail region is noticeable. On the other hand, the prediction by 
the untruncated normal model (dashed line) differs substantially from 
the data. The equations and the calculation procedure of the gamma 
model have been described in Ref. 7. The estimated statistical param- 
eters of apparatus-case crosstalk based on the simple gamma model 
are listed in Table I. 

Many authors have used the Wilkinson’s method” to calculate the 
power sum distribution from the pair-to-pair distribution. The obvious 
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Table |—Statistical parameters of apparatus-case 
crosstalk 
(Eight 466-type apparatus cases measured in Atlanta) 


Gamma Model} (For all data) 


yx 3.30 x 107-8 

Oy 2.23 x 10°° 
Pair-to-pair us B00 

B 0.833 

M,(dB) 96.0* 

o, (dB) 10.7* 

& 1.65 x 1077 

Os 1.58 x 1077 
Power sum of 50 pair-to-pair mM 568.00 

crosstalk losses a 6.58 
Moa(dB) 69.30 
0g(dB) 3.62 


* These values are estimated by gamma model and are slightly different 
from those of normal model. 

+ The definitions of terms and equations related to gamma model are given 
in Ref. 7. 


discrepancy between the measured data and the dashed line predicted 
by the untruncated normal model in Fig. 5 has prompted some authors 
to abandon the Wilkinson’s method” entirely and to simply fit the 
measured power sum distribution on Figs. 3 and 5 by a normal 
distribution. As will be shown later, this approach is too optimistic by 
4 dB at the critical 0.1 percent point. Thus, the normal model faces a 
dilemma of being too pessimistic (see Fig. 5), if Wilkinson’s method of 
power sum calculation is used, and being too optimistic at the 0.1 
percent point, if Wilkinson’s method is by-passed (i.e., simply fit the 
measured power sum distribution by a normal distribution). The use 
of truncated normal model for pair-to-pair distribution suffers a draw- 
back of uncertain truncation point as discussed in Ref. 7 and several 
dBs of error at the 0.1 percent point just mentioned. 

In engineering applications, the behavior of the power sum distri- 
bution in the lower tail region is most important because the engi- 
neering objective of T1 systems is set at the 0.1 percent point for 50- 
section metropolitan applications. Unfortunately, Figs. 3 and 5 show 
that the measured data deviate substantially from the predictions by 
gamma and normal models in the important lower tail region. These 
discrepancies are predictable by both the Monte Carlo and the lower 
bound method as discussed in the next section. 

Reference 7 and this paper indicate that an accurate prediction of 
crosstalk power sum distribution from pair-to-pair distribution is often 
difficult. One is tempted to abandon the pair-to-pair crosstalk statistics 
entirely and to rely solely on the measured power sum distribution. 
However, the studies of pair-to-pair statistics and other decomposi- 
tions, such as within-slot versus non-within-slot crosstalk, and within- 
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harness versus non-within-harness crosstalk, are necessary to provide 
insights for understanding and for proper modeling of non-Gaussian 
power sum distribution. The technique for predicting power sum 
distribution from pair-to-pair distribution is also necessary in charac- 
terizing some practical situations where the apparatus cases are only 
partially filled. 


lll. POWER SUM CALCULATIONS BY MONTE CARLO AND LOWER 
BOUND METHODS 


The extensive laboratory and field measurements indicate that the 
distribution of the power sum of apparatus-case crosstalk has a notice- 
able bend towards more severe crosstalk levels in the lower tail region 
as shown in Fig. 3. The lower tail region of the measured data on Fig. 
3 has an effective standard deviation (i.e., slope) of 6 dB on the normal 
probability coordinates. This slope agrees very well with the slope of 
the measured distribution of T1 repeater section margins in the lower 
tail region.” Thus, the laboratory measurements and field measure- 
ments consistently indicate that the distribution of the power sum of 
ACXT cannot be described precisely by a simple model, such as the 
normal or the gamma distributions. 

With such understanding, we will avoid assuming any simple model 
for the total power sum distribution and will use more sophisticated 
techniques, such as the Monte Carlo method or the lower bound 
technique to obtain the correct power sum distribution. 

Previous studies by Marlow’ and Janos” indicate that the power 
sum distributions will have a bend if there are strong, dominant 
components whose mean or standard deviation differs substantially 
from those of the other components of the power sum. Under such 
circumstances, the lower tail of the power sum distribution will behave 
like that of the dominant components and, hence, show a bend. With 
this hint, we naturally look for the possible existence of dominant 
components in the power sum of apparatus-case crosstalk disturbers. 

Figures 1 and 2 show the repeater slot and wiring arrangements of 
T1 apparatus case. Each T1 system in an apparatus case suffers from 
two within-slot disturbers* and 48 non-within-slot disturbers assuming 
the case is 100 percent filled. The Atlanta laboratory data show that 
the mean value of the within-slot crosstalk is 15 dB worse than that of 
the non-within-slot crosstalk. Such a large difference means that the 
within-slot crosstalk and the non-within-slot crosstalk must be treated 
separately in the power sum calculations. 

The circles on Fig. 6 show the distribution of power sum of the 48 
non-within-slot crosstalk disturbers measured in Atlanta Laboratory. 


*Each apparatus case slot holds two T1 regenerators. 
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Fig. 6—Distribution of power sum of 48, pair-to-pair, non-within-slot, crosstalk losses 
of 466-type Tl apparatus case. Data measured from eight apparatus cases in Bell 
Laboratories, Atlanta. Data represents between-slot crosstalk component of data in 
Fig. 5. 


The circles, triangles, and crosses on Fig. 7 show the distributions of 
power sum of the two within-slot-crosstalk disturbers measured in 
Atlanta, Illinois, and California, respectively. The solid lines on Figs. 
6 and 7 are the corresponding gamma approximations with the param- 
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Fig. 7—Distributions of power sum of two, pair-to-pair, within-slot, crosstalk losses of 
466-type T1 apparatus case from laboratory measurements in Atlanta and field mea- 
surements in Illinois and California. Data represents within-slot crosstalk component of 
data in Figs. 3 and 5. 
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Table I[—Statistical parameters 
of non-within-slot crosstalk of 


apparatus case 
(Eight 466-type apparatus cases 
measured in Atlanta) 


Gamma model for power sum of 48 
non-within-slot crosstalk (Fig. 6) 





5 8.99 x 10°° 
Os 9.10 x 10°8 
im 551.00 
a 6.20 
Mo(dB) 72.20 
0g (dB) 3.79 


eters listed in Tables II and III, respectively. The total power sum of 
the 50 pair-to-pair crosstalk disturbers is equal to the power sum of 
the two gamma distributed random variables on Figs. 6 and 7. Tables 
II and III show that the mean values of these two gamma variates 
differ by only 7 percent (i.e., 72.2 vs. 77.5 dB), whereas the standard 
deviations differ by 100 percent (i.e., 3.8 versus 8.0 dB). A bend on the 
distribution of their power sum is, therefore, expected. 

The Monte Carlo method for power sum calculation has been used 
by many authors.'”"® The lower bound technique is described in the 
Appendix. Figure 3 shows that both the Monte Carlo and the lower 
bound methods predict the bend in the total power sum distribution 
in close agreement with the measured data. The Monte Carlo result* 
agrees very well with the measured data over the entire range. The 
predicted lower bound (in probability) is practically identical to the 
Monte Carlo result in the lower tail region (<2 percent) and is 
applicable to T1 system engineering. The lower bound method has the 
advantages of being simple computationally and of providing some 
physical insights into the power sum behavior in the tail region as 
discussed below. 

Let x denote the power sum resulting from within-slot crosstalk 
(Fig. 7) and let y denote the power sum due to non-within-slot crosstalk 
(Fig. 6). Furthermore, let z denote the power sum of x and y, the total 
power sum of 50 pair-to-pair crosstalk disturbers. The Appendix shows 
that a lower bound, Pza(z S 5), of the probability distribution of z is: 


Prz(z S 6) = P(x Sb) + P(y sb) —P(x=b)-P(ysb), (1) 


where 6 represents the crosstalk level at which the probabilities are of 
interest. The data in Figs. 6 and 7 show that 


P(iysb)<P(x=b) for b= 62 dB, (2) 


*A sample size of 5000 is used in obtaining the Monte Carlo result (the solid line) in 
Fig. 3. 
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Table IIl—Statistical 
parameters of within-slot 


crosstalk of apparatus case 
(466-Type apparatus cases measured 
in Atlanta, Illinois, and California) 


Gamma model for power sum of two 
within-slot crosstalk (Fig. 7). 


s 8.45 x 107° 
Os 2.75 x 107? 
p 101.00 

a 1.258 
Ma(dB) 77.50 

og (dB) 8.00* 


* Since a gamma distribution is not a 
straight line on a normal probability co- 
ordinates, the slope of the gamma distri- 
bution in the lower tail region is not the 
same as that (i.e., 8 dB of o) in the middle 
range. 


which implies that: 
Pria(z sb) = P(x <b) for b= 62GB. (3) 


The data in Figs. 3 and 7 indeed support this simple approximation. 
Therefore, the lower bound method demonstrates through eq. (3) that 
the total power sum distribution behaves like that of the dominant 
component x in the lower tail region (i.e., for those situations where 
ACXT is worse than 62 dB). Notice that the dominant component x 
represents the power sum of the two within-slot-crosstalk losses. 


Vi. CONCLUSION 


The extensive data on T1 apparatus-case crosstalk for 466-type 
cases from laboratory measurement in Atlanta and field measurements 
in Illinois, California, and Texas consistently indicate that the power 
sum distribution has a bend towards more severe crosstalk levels in 
the lower tail region. It is shown that this bend is because of the 
dominant effect of within-slot crosstalk. Both the Monte Carlo and 
the lower bound methods of power sum calculations predict this bend 
if the power sum contains a dominant component which differs sub- 
stantially from other components. 

The 0.1 percentile of the distribution of power sum of ACxT is about 
55 dB for 466-type apparatus case. This is about 4 dB worse than that 
predicted by the conventional normal model which ignores the bend 
in the tail region. 
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APPENDIX 
A Lower Bound of Power Sum Distribution 


This appendix describes a technique to obtain a lower bound of the 
probability distribution of crosstalk power sum. In the low probability 
tail region, this lower bound is fairly tight and provides a simple 
approximation to the power sum distribution. This approach is inspired 
by the work of Marlow and Farley.’®”° 

Let P(x = b) and P(y S bd) be the cumulative distributions of two, 
positive, independent, random variables x and y, respectively, and let 


z=-10 logif 10 + 10} (4) 
be the power sum of x and y. This definition implies that 
z= min(x, y), (5) 
and 
P(z = b) = P(min(x, y) = b) (6) 
=P(x=b,y=b), © (7) 


where min(x, y) denotes the minimum of x and y, and P(x = b, y = b) 
denotes the probability that both x and y exceed b. The independence 
of x and y implies that 


P(x = b, y= b) = P(x = b)-P(y= B). (8) 
By definition: 
P(z< b) =1-P(z = b) 
P(x < b) =1— P(x 2 b) 
P(iy<b)=1-P(y=b) . (9) 
Combining eqs. (7), (8), and (9) yields 
P(z < 6) =1—- P(x= b)-P(y= d) 
=1-[1-—- P(x < b)]-[1 — P(y < b)] 
= P(x < 6b) + P(y <b) — P(x < b)-P(y< 5b). (10) 
Therefore, the right-hand side of eq. (10) represents a lower bound (in 
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probability) for the distribution of the power sum z. Notice that this 
lower bound can easily be calculated when the distributions P(x < b) 
and P( y = b) of the components x and y are known. 

In the lower tail region where both P(x = b) and P(y = db) are small, 
the probability of both x and y being less than 6b simultaneously is 
extremely small. Therefore, the inequality eqs. (5) and (10) asymtoti- 
cally approach equalities in the lower tail region (i.e., when 6 is small). 
This means that in the low probability tail region, the lower bound eq. 
(10) is fairly tight and provides a simple but accurate approximation 
to the power sum distribution. 
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Current-Carrying Capacity of Fine-Line Printed 
Conductors 


By A. J. RAINAL 
(Manuscript received February 23, 1981) 


This paper presents simple equations, along with experimentally 
determined parameters, to calculate the transient temperature rise of 
current-carrying, fine-line (~7 mils) conductors on various styles of 
circuit packs. The styles of circuit packs include wire wrap, double- 
sided, metal, bonded, and various multilayer boards. All styles of 
circuit packs in the BELLPAC™ system family are included. The 
maximum steady-state temperature rise of a nicked or constricted 
current-carrying conductor is also treated. The calculated transient 
and steady-state temperature rises agree with experimental results. 


|. INTRODUCTION 


Printed wiring technology presently provides the physical designer 
with fine-line copper conductors to interconnect integrated circuits 
and other components at the circuit-pack (cp) level. These fine-line 
conductors now have a nominal width of 7 mil, and a nominal thickness 
of 1.4 mil. For such relatively small conductor sizes, (~AWG39), the 
current-carrying capacity of the conductors becomes an important 
matter of concern. Can such fine-line conductors carry the required 
current to operate the various components on a CP without causing an 
excessive temperature rise? What is the temperature rise during nor- 
mal current flow? If a fault occurs and a current of 10 A flows for 100 
ms, will the cP be damaged? Questions of this nature are becoming 
very important as printed wiring technology provides finer conductors 
for the electrical interconnections. 

Also, as large-scale integrated (LSI) circuits are introduced, the 
assembled cps are becoming more expensive. Therefore, there is more 
incentive to protect the cp from possible damage as a result of an over 
current. 

Some early work by W. Aung and A. J. Colucci’ has shown that a 
7.5-A current flowing through a fine-line printed conductor inside a 
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multilayer board (MLB) can cause the MLB to break into flames in 
about 5 s. Also, W. T. Smith” reported that a current flow through a 
fine-line printed conductor (2.8-mils thick) should be limited to about 
5 A if temperature rises are to be limited to 120°C. 

Clearly, to avoid possible disaster or damage, the physical designer 
of electronic equipment must be able to estimate the transient tem- 
perature rise of fine-line, current-carrying conductors on all styles of 
circuit packs and backplanes. 

During the past few years, the Bell System has introduced a modular 
packaging system (BELLPAC* system’) for packaging electronic 
equipment. This system makes use of a number of cp styles that have 
common features suitable for computer-aided design. Fine-line printed 
conductors are available for use on any of the cP styles. The BELLPAC 
system project has provided us with the opportunity to study the 
current-carrying capacity of fine-line conductors on a variety of CP 
styles. 

The purpose of this paper is to present some useful results concern- 
ing the current carrying capacity of fine-line conductors on various 
styles of cps. The results include the transient temperature rise of a 
current-carrying printed conductor, and the maximum steady-state 
temperature rise of a conductor which may be nicked or constricted. 

For the special case of a double-sided epoxy printed wiring board 
(flex or rigid), some results concerning the steady-state temperature 
rise of fine-line printed conductors have been reported in Refs. 4 and 
5. Although the methods reported in these two references differ, the 
results agree well with one another. 

A listing of the cp styles of interest in this paper, along with a short 
description of each is presented in Table I. Figure 1 shows the corre- 
sponding physical layups of the cp styles. These physical layups 
include all of the cp styles presently in the BELLPAC system. We 
shall see that the results in this paper can be applied to estimate the 
current-carrying capacity of fine-line conductors on any layer of any of 
the cp styles shown in Fig. 1. 


il. BASIC EQUATIONS 
2.1 General results 


For a general current-carrying conductor, the conservation of heat 
energy requires that the average temperature rise satisfy the following 
differential equation: 


ee — AT 
PRi{1 + a AT Jat = CdaT + 5 dt, (1) 
T 


* Trademark of Western Electric. 
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Table I—Description of the circuit-pack styles 


Circuit-Pack Style Description 
Wire wrap Wire-wrap board for breadboarding 
Extender board 6-Layer MLB, 2-pad layers, 2 signal layers, power (P) 


and ground (G) on inside, dedicated ground conduc- 
tor between every pair of signal conductors 


Double-sided (epoxy) Double-sided, epoxy PWB 
Double-sided (metal) Double-sided, metal core, PWB 
Bonded board (LAM- Flex bonded to epoxy-coated steel 
PAC)* 
4L MLB (EXT P/G) 4-Layer MLB, 2 signal layers, P and G on outside 
6L MLB (EXT P/G) 6-Layer MLB, 4 signal layers, P and G on outside 
6L MLB (INT P/G) 6-Layer MLB, 2 pad layers, 2 signal layers, P and G on 
inside 
6L MLB (INT P/G, Sur- 6-Layer MLB, 4 signal layers, P and G on inside 
face Routing) 
8L MLB (INT P/G) 8-Layer MLB, 2 pad layers, 4 signal layers, P and G on 
inside 


* This particular bonded board is also known as LAMPAC. 


where 


I = current flow through the conductor (amperes) 

R, = resistance of the conductor at ambient temperature (ohms) 

a, = ambient temperature coefficient (per degree C) 

= [T, + 234.45]"' (a good approximation in the case of copper) 
T, = ambient temperature (°C) 
t = time duration of the current flow (seconds) 

C = thermal capacity of the material heated (J/°C) 
Rr = thermal resistance of the conductor on the cp style of interest 
__ (°C/W) 
AT = average temperature rise (°C) along the length of the conductor. 


The left-hand side of eq. (1) represents the energy dissipated, and the 
two terms on the right-hand side represent the stored and radiated 
energy, respectively. More general forms of eq. (1) are discussed in 
Ref. 6. 

Laboratory measurements show that the thermal capacity, C, is 
time dependent. The physical reason for this dependence is that more 
and more of the material is heated as time goes on. Initially, only the 
conductor and a small portion of the substrate and covercoat are 
heated. 

However, if the time axis is partitioned into appropriate time inter- 
vals, it turns out that the thermal capacity, C, is approximately 
constant over each of the time intervals. The solution to eq. (1) in the 
case of three such time intervals is given by: 


AT = ATs (1 — exp{—(RrQi) (1 —?RiRray)t}] OSt<th (2) 
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Fig. 1—Lay-ups of the various circuit-pack styles. 


1378 THE BELL SYSTEM TECHNICAL JOURNAL, SEPTEMBER 1981 


AT = AT t= exp| - —(RrC,) (1 — P-RiRren) 


Cu Ci 
. ~_-— t hZstst. 
(: C, Cy =) 4] ae (3) 


AT = ATs E - exp] ~(rC)¥0 — I?RiRra;) 


Oy Cit, Cit At 
-— es te<t< 0, 
(1-2 Ce Cn GC: *) | : = 








(4) 
where 


— PRiRr : 
AT,; => =~ = Steady-state average temperature rise, 
(1 — P-RiRras) : aan 
C, = thermal capacity during the time interval [0, f:], 
C2 = thermal capacity during the time interval [t, ¢2], 


C3 = thermal capacity during the time interval [f2, ©]. 


In general, one can partition the time axis into n contiguous time 
intervals and obtain a set of n equations. For our purposes, n = 3 
proved to be sufficient. 

The current-carrying capacity of a conductor is limited by the 
permissible temperature rise of the conductor above the ambient 
temperature. Therefore, once the pertinent parameters are known, the 
above equations can be used to calculate the current-carrying capacity 
of a particular conductor. 

In Section ITI, we describe the experimental method used to measure 
all of the pertinent parameters. 


2.2 Some special cases 
2.2.1 Runaway or critical current 


The functional form of AT, shows that a runaway or critical current, 
I,, exists for which AT,, — 0. That is, as I > I., AT,, —> ©, and the 
current-carrying conductor never reaches steady-state temperature. 
The value of J, is given by: 


baa (5) 
(Rate ———— 

Rikray 
The phenomenon of runaway and the value of runaway current is 
consistent with our experience in the laboratory. We found that as we 
approached the critical value, [., the temperature of the conductor 
rises rapidly beyond the tolerable limits of the substrate and perma- 
nent damage results. 


2.2.2 Small t — initial temperature rise 
From eq. (2), we find that as t > 0, we have 


FINE-LINE PRINTED CONDUCTORS 1379 


set ER yt 
AT-= ‘ 
Ci 


Notice that this result is independent of the thermal resistance Rr. 
That is, the initial heating process is adiabatic. 





(6) 


2.2.3 Large t — steady state temperature rise 
From eq. (4), with I < I., we see that as t > o we have 
PRiRr 


AT HAT = 
(1 = IP?Ri\Rra) 


(7) 

Equation (7) shows that the steady-state temperature rise depends 
on the product of the electrical resistance (a property of the conductor) 
and the thermal resistance (a property of the environment of the 
conductor). For a given current J, and ambient T:, one can only reduce 
AT:;; by reducing the product RiRr. 


lll. EXPERIMENTAL DETERMINATION OF THE PARAMETERS 


In order to carry out this study, appropriate test boards were 
designed for each cp style shown in Fig. 1. Except for the double-sided 
metal board, all test boards were fabricated at the Western Electric 
printed-circuit manufacturing plant at Richmond, Virginia. The dou- 
ble-sided metal board was manufactured at the Western Electric plant 
in Kearny, New Jersey. 

In all cases, the physical dimensions of the printed conductors had 
the nominal values of length L = 12 in., width W = 7 mil, and copper 
thickness f = 0.5 to 3 mil. 

The experimental method used to determine the pertinent param- 
eters is based on measuring, indirectly, the average temperature rise, 
AT, along the conductor as a function of time when a step function of 
current is applied. The measurement is based on the well-known 
resistance thermometer formula:®*’ 


_ R= Rl + wAT, (8) 
where 


V = voltage across the conductor, 
I = magnitude of the step function of current, 
R = measured resistance of the conductor. 


The procedure used to measure AT is as follows: The ambient 
temperature JT is recorded and A; is measured by means of an ac 
Kelvin bridge. This measurement involves a small current (100 mA or 
less) which causes a negligible temperature rise. Then a step function 
of J amperes is directed through the conductor of interest, and the 
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resulting voltage drop, V, across the conductor is recorded as a function 
of time. Equation (8) is then used to deduce the corresponding AT' as 
a function of time. 

To help ease the data gathering, a Kaye Instruments digistrip 
transmitter was used to format the data for print out and magnetic 
tape storage on a Texas Instruments model 733 data terminal. Subse- 
quently, the data was transmitted (via an acoustical coupler) over the 
telephone line to the computation center for storage. At this point, 
special computer programs were used to edit the stored data and to 
produce the computer plots. 

The procedure used to determine the constants appearing in eqs. 
(2), (3), and (4) is as follows: The average temperature rise is first 
computed from eq. (8) by using the steady state voltage value corre- 
sponding to a current flow of J amperes for a sufficiently long time 
(usually about 10 min). This is repeated for a number of different 
values of current. Then, the best value (minimum mean-square-error 
sense) of the product RR, is determined from the measured data and 
eq. (7) which can be rewritten as 

AT ss = R, RI’. (9) 
iL + aATss 





If the left-hand side of eq. (9) is plotted as a function of J’, the slope, 
Rif, of the best fitting line is the quantity of interest. Since the value 
of R, is known, Rr can be determined. 
According to eqs. (2), (3), and (4), a plot of the measured values of 
Y= -In| 1 a ad 
AT;s 





versus time, ¢, yields points which tend to fit a series of approximately 
broken lines of positive slope. From this plot, values of f; and f2 can be 
selected as the break points of these broken lines. In the region 0 = ¢ 
= t,, the best value of Ci (minimum mean-square-error sense) is 
determined by equating the slope of the best fitting line to the slope of 
the negative of the exponent in eq. (2). In a similar manner, C2 is 
determined by using the measured data and the slope of the negative 
of the exponent in eq. (3). Finally, C3 is determined by using the 
measured data and the slope of the negative of the exponent in eq. (4). 
At this point, all of the parameters needed in eqs. (2), (3), and (4) are 
known, and these equations can be used to calculate the average 
temperature along the conductor as a function of applied current, time, 
conductor resistance, and ambient temperature. 

In this manner, the pertinent parameters were determined for fine- 
line conductors on all signal layers of all cp styles shown in Fig. 1. The 
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resulting parameters are presented in Table II. The values of ¢#; = 0.55 
s and f2 = 3.55 s were found to apply to all of the cp styles. The results 
were scaled to a conductor length of 12 in. and width of 7 mil by using 
the scaling laws C; ~ L, Rr ~ 1/L, and H ~ 1/W. 

The values of H listed in Fig. 2 and Table II are steady-state 
parameters and will be discussed in more detail in Section V. 


IV. EXPERIMENTAL VERIFICATION 


Figure 2 presents the experimental values of Y(t) for the case of a 
double-sided epoxy printed wiring board (PWB) with covercoat. The 
parameters t), tz, Ci, C2, Cs, and Rr were determined by the methods 
described in Section III. The final C; were determined by averaging 
the results over three different values of current. 

For the double-sided epoxy cp, Fig. 3 compares the experimental 
values of the transient temperature rise, AT(t), with the corresponding 
theoretical values. The theoretical values were determined by using 


2.5 


t7 =0.55 s 

t2 =3.55 s DOUBLE-SIDED (EPOXY) 
°C LAYER G/So 

Rr= 13.94 Ww WITH COVERCOAT 


Ww 
2.0. 4 = 0.506E — 06 


J 
C7 = 0.093 — 


J 
C2 = 0.61 2 


J 
— C3 =3.27 — 
1.5 3 ae 





L = 11.80 in 
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" 

ss 
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oO /=2.75 A (ATgs = 68.2 °C) 
4 1=3.00 A (ATs = 81.0 °C) 
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Fig. 2—Experimental values of Y(t). Slopes of the broken lines determine the average 
thermal capacities, Cj, in the three time intervals. 
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Table II—Measured parameters for various circuit-pack styles 
(Conductor length = 1 ft, conductor width = 7 mil, and wire-wrap 
diameter = 10 mil) 


J J J °C Watts 
Circuit-Pack Style* Owq Gag Gag Arp mil? °C 
Wire wrap 
(Milene)+ 0.118 0.170 0.265 30.5 0.869 1077 
(Teflon)+ 0.150 0.259 0.406 30.2 0.877 1077 
Extender board 0.121 0.697 8.71 6.32 0.942 10° 
Double-sided (epoxy) 0.095 0.623 3.32 13.71 0.434 107° 
Double-sided (metal) 0.052 0.3805 11.22 9.00 0.661 10~* 
Bonded board P/S1 0.083 0.311 7.26 8.27 0.720 10°° 
(LAMPAC) G/S1 0.090 0.415 9.07 7.35 0.810 10-° 
4L MLB (ext P/G) 0.124 1.06 3.65 11.49 0.518 10°® 
6L MLB (ExT P/G) 
Si, S4 0.124 1.06 3.65 11.49 0.518 107° 
0.179 1.38 4.75 9.55 0.623 10°° 


So, S3 
6L MLB (Int P/G) 


ly D2 0.121 0.697 8.71 6.32 0.942 10°° 
6L MLB (nt P/G, Surface 


Routing) 
Si, S4 0.076 0.443 5.59 10.47 0.569 10°® 
So, S3 0.112 0.887 7.61 7.33 0.812 10~¢ 
8L MLB (nT P/G) 
Si, S4 0.076 0.438 5.01 9.70 0.614 10-® 
So, S3 0.108 0.888 6.87 7.38 0.806 107° 


* The circuit-pack styles having surface conductors were covercoated. 
Trademark of W. L. Gore & Associates, Inc. 
Trademark of E. I. DuPont de Nemours, Inc. 


the constants 4), t2, Ci, Co, Cs, Rr in eqs. (2), (3), and (4). Figure 3, 
shows that this method of estimating the transient temperature rise 
agrees with experimental results. 

Similar plots have verified that this method of estimating the 
transient temperature rise also agrees with experimental results for all 
other cases of interest in this paper. 


V. SOME APPLICATIONS 
5.1 Steady state temperature rises 


The steady-state temperature rises of fine-line printed conductors, 
or wire-wrap conductors, can be readily computed from eq. (7). The 
required values of thermal resistance, R7z, are listed in Table II for a 
conductor length of 1 ft, printed conductor width of 7 mil, and wire- 
wrap diameter of 10 mil. Also, the required values of the electrical 
resistance, #, can be computed from: 


pL . 
RS 
1 TW (printed conductor), (10) 
or 
4oL 
Rk, = i (wire-wrap conductor), (11) 
77. 
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EXPERIMENTAL VALUES 

O / = 2.50 A (ATs = 55.1 °C) 
/ = 2.75 A (ATs = 68.2 °C) 
1 = 3,00 A (A755 = 81.0 °C) 
THEORETICAL VALUES 


THEORETICAL AT, = 83.90 °C 
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—= 
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Fig. 3—Experimental values of the transient temperature rise, AT(t), are compared 
with theoretical values. 


where 


p = (0.67878) 10~*[1 + 0.00393(7; — 20)] ohm-mil 
T, = ambient temperature, °C 

L = length of conductor = 1 ft (12,000 mil) 

t) = thickness of conductor, mil 
W = width of conductor = 7 mil 

D = diameter of wire-wrap conductor = 10 mil. 


As an example, consider a double-sided (epoxy) cP style. From Table 
II, we see that Rr = 13.71°C/W. From eq. (10), for L = 12,000 mils, 
t) = 1.4 mils (1 oz cu), W = 7 mil and 7; = 20°C, we find that RR; = 
0.8312. For a current flow of J = 2.5 A, eq. (7) then yields AT,, = 
98.9°C. 

For conductor lengths other than L = 1 ft, one can use the scaling 
law Rr ~ 1/L to scale the values of Rr listed in Table II. Also, Rr is 
essentially independent of W as is shown by eq. (16) of Ref. 4. 

Let us now compute the maximum steady-state temperature rise, 
max AT, when the same current-carrying conductor is nicked or 
constricted. The maximum temperature rise of such conductors was 
treated in Ref. 4. It was shown that the key parameter was the value 
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of H, the coefficient of surface heat transfer. The values of H for the 
cP styles of interest in this paper are listed in Table II. 

As an example of the effect of a nicked or constricted conductor on 
steady-state temperature rise, Fig. 4 compares the computed results 
for a fine-line current-carrying conductor on a double-sided (epoxy) 
board and a 6-layer MLB (INT P/G, surface routing). At a current of 
I = 2.5 A, the nicked conductor on the double-sided board rises in 
temperature to about 119°C; whereas, the nicked conductor on the 
inside signal layer (S_) of the 6-layer MLB rises in temperature to only 
56°C. 

Tabulated results for the special case H = (0.52)107° (Watts/mil?°C) 
were presented in Ref. 4. Since many of the values of H tabulated in 
Table II are close to this value, the earlier results can also be applied 
to many of the cp styles of interest in this paper. Also the simple 
relationship given as eq. (13) of Ref. 4 yields results which agree well 
with those presented in Fig. 4. 

Table II gives the values of H when W = 7 mil. For conductor 
widths other than W = 7 mil, one can use the scaling law H ~ 1/W. 


DOUBLE-SIDED EPOXY /L¢ = 2W 
W = 7mil Lo =0 
We = 3.5 mil 
L = 12000 mil 
Le = 0,14 mil 
Ty = 20°C 
to =1.4 mil 


MAXIMUM AT¢ (°C) 





0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 
7 IN AMPERES 


Fig. 4—The effect of a nicked or constricted conductor on the maximum steady state 
temperature rise. 
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5.2 Transient temperature rises 


The transient temperature rises of fine-line printed conductors, or 
wire-wrap conductors can be computed from eqs. (2), (3), and (4). The 
required values of C; and Rr are listed in Table II, also the values of fy, 
and f, are given by «4 = 0.55 s, and #2 = 3.55 s, as was discussed in 
Section ITI. 

An example of the transient temperature rise of a fine-line printed 
conductor on a double-sided cp is presented in Fig. 5. For this case, 
the appropriate parameters are listed in Table II as: 


°C 
Watt 





J 
C2 = 0.623 a C3 = 3.32 —, and Rr = 13.71 


J 
C; = 0.095 = oC oG 


it (Bd 
5.3 Temperature rises resulting from fault currents 


A signal conductor on a cP normally carries a maximum current of 
about 0.1 A. Figure 4 shows that the temperature rise of such a current- 
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Fig. 5—Theoretical values of the transient temperature rise, AT \(t), are compared for 
different styles of circuit packs. 
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carrying conductor is negligible. However, a component failure or 
other malfunction can result in a current flow of many more amperes 
through the fine-line signal conductor. In this case, a circuit fuse or 
other over-current protection device can take many milliseconds to 
interrupt the flow of current. In such situations, the results in this 
paper can be applied to calculate the temperature rises on any layer of 
any of the cp styles shown in Fig. 1. 

For example, consider the case of a signal conductor having the 
physical dimensions of L = 1 ft, W = 7 mil, and f& = 1.4 mil ona 
double-sided (epoxy) cp. Let us also assume that the ambient temper- 
ature is 50°C. Suppose the fault current is 10 A and the fuse or over- 
current protection device interrupts the fault current in 100 ms. What 
is the resulting temperature of the current-carrying conductor imme- 
diately before the current flow is interrupted? 

From eq. (10), the electrical resistance of the conductor is R; = 
0.9292. From eq. (6) we find that 


a 2 2 
apn PR _ (10. A%09290)018) _ grag ay 


The value of C; was taken from Table II. Thus, the temperature of the 
signal conductor is about 50°C + 98°C = 148°C. From some additional 
experimental work, we have found that the epoxy glass substrate 
begins to discolor at about 175°C. Thus, in our example, the substrate 
would not be discolored if the fault current of 10 A is interrupted in 
about 100 ms. Of course, in an application, one may want to restrict 
the conductor temperature to much less than 175°C. 

It is important to notice that the same conductor temperature would 
result even if the length, L, of the conductor were only a few inches, 
since both R, and C; are proportional to L. Also, for conductor widths 
other than W = 7 mil, one can use the scaling law C,; ~ W to scale the 
values of C, listed in Table IT. 


5.4 Some comparative results 


Figure 5 compares the transient temperature rises of fine-line printed 
conductors carrying a current of J = 2.5 A on various cP styles. In the 
steady state, the wire-wrap conductors exhibit the least rise in tem- 
perature (21°C); whereas, the conductor on the double-sided cp yields 
the highest rise in temperature (99°C). 

Notice that the relative current-carrying capacity of the various CP 
styles depends on the duration of current flow. For example, at about 
t = 5 s, the conductor on the double-sided metal board exhibits the 
highest temperature rise (47°C). Interestingly, this result shows that 
the thermal conductivity of the dielectric resin of the metal board is 
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somewhat less than the thermal conductivity of epoxy glass. Similarly, 
Fig. 5 also shows that Milene insulation has a somewhat lower thermal 
conductivity than Teflon insulation. 

In most applications, one is usually interested in the comparative 
results during the steady state. In this case, the double-sided (epoxy) 
board is the worst from the point of view of current-carrying capacity. 

Finally, notice that in all cases, the conductor temperature rises to 
a substantial fraction of its final value in the first five seconds. 


Vi. SUMMARY 


This paper presents simple equations, based on the conservation of 
heat energy, to calculate the transient temperature rise of current- 
carrying, fine-line (~7 mil) conductors on various styles of circuit 
packs. The equations depend on various parameters which were de- 
termined experimentally. All of the styles of circuit packs in the 
BELLPAC system are included. These styles range from wire wrap to 
various MLBS. 

The maximum steady-state temperature rise of a nicked or con- 
stricted current-carrying printed conductor is also treated. 
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A Comparative Study of Several Dynamic 
Time-Warping Algorithms for Connected-Word 
Recognition 


By C. S. MYERS and L. R. RABINER 
(Manuscript received February 11, 1981) 


Several different algorithms have been proposed for time register- 
ing a test pattern and a concatenated (isolated word) sequence of 
reference patterns for automatic connected-word recognition. These 
algorithms include the two-level, dynamic programming algorithm, 
the sampling approach and the level-building approach. In this 
paper, we discuss the theoretical differences and similarities among 
the various algorithms. An experimental comparison of these algo- 
rithms for a connected-digit recognition task is also given. The 
comparison shows that for typical applications, the level-building 
algorithm performs better than either the two-level DPp-matching or 
the sampling algorithm. 


1. INTRODUCTION 


Research in the area of automatic speech recognition has progressed 
to the point at which a wide variety of isolated word recognition 
systems have been implemented and used successfully for many ap- 
plications."? These systems have been used in such applications as 
data entry, searching, and sorting; however, use of these systems is 
restricted by the format of the speech input, i.e. isolated words. For 
many applications, a connected-word input format would have several 
advantages. Examples of such applications include: 

(t) Credit card entry—A sequence of digits (and possibly letters) 
is required to specify the credit card number. 
(tt) Directory listing retrieval—A sequence of letters is used to 
spell the name for which a directory listing is required. 

(zit) Airline reservations—Sentences based on a restricted vocabu- 

lary and a restricted syntax are used to make reservations. 
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Several techniques for recognizing connected-word sequences from 
isolated word reference patterns have recently been proposed.’*” 

In this paper, we compare, both theoretically and experimentally, 
three algorithms which have been proposed for connected-word rec- 
ognition. These algorithms are the two-level, dynamic programming 
matching (TLDPM) approach, the sampling approach, and the level- 
building (LB) approach.’*” Another algorithm of the same general 
class as the ones considered here has been recently proposed by Bridle; 
however, it is not considered here.’* 

It should be noted that all of the algorithms which have been 
proposed for connected-word recognition are loosely related to a 
general information-theory-based algorithm proposed by Bahl and 
Jelinek.” However, each of these connected-word recognition algo- 
rithms make different assumptions, have different implementations, 
and make specific tradeoffs. Thus, the properties of these algorithms 
are entirely different from those of Bahl and Jelinek; therefore, they 
warrant independent study. 

In Section IJ, we review each of the three connected-word recogni- 
tion algorithms and then provide a theoretical comparison of the 
general properties. In Section III, we theoretically compare the con- 
nected-word recognition algorithms, in Section IV, we describe and 
give results of several experimental comparisons of these algorithms, 
and in Section V, we discuss these results and their implications for 
practical word-recognition systems. 


ll. CONNECTED-WORD RECOGNITION ALGORITHM 

Figure 1 shows the basic structure for all the connected-word rec- 
ognition algorithms under consideration. The feature extraction is 
generally similar to that used in most isolated word-recognition sys- 
tems. Typical feature sets include energy of a set of bandpass filters,”® 
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Fig. 1—Block diagram of a generic connected word recognition system. 
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and LPC coefficients.’* Following feature extraction, the test utterance, 
now represented by a sequence of frames of the feature vector, is 
nonlinearly time-registered, with a set of reference patterns, and a set 
of distance, or dissimilarity, scores are calculated. Following time 
registration, a decision rule chooses the sequence of reference patterns 
which best matches the test utterance. Feedback, in the form of 
positional information, is used to determine which portion of the test 
utterance is to be matched to a given reference pattern. It is generally 
assumed that the test utterance consists of a sequence of words spoken 
by a cooperative user; that is, the words of the string are carefully 
articulated and spoken in a slow, deliberate manner, but not as isolated 
words. However, the reference patterns do consist of isolated words. 
The goal of the connected-word recognizer is to find the sequence of 
concatenated reference patterns, subject to given syntactical con- 
straints, that best matches the test pattern. Among the issues which 
arise in solving for the best string are the following: 

(tz) How can the reference and test patterns be time-registered? 

(tt) How can the reference patterns be modified to account for 
both coarticulation and the natural shortening of words inherent in 
connected speech? 

(wi) Is the determination of the best concatenation of reference 
patterns done in a sequential manner, i.e. one decision at a time, or is 
some form of backtracking allowed? 

(tv) How can alternative matches to the test utterance be generated 
in addition to the best match? 

(v) Can syntactical constraints be used explicitly in the recognition 
stage or is some form of post-processing required? 

In the following, we review the three connected-word recognition 
algorithms and discuss how each of them answers these questions. 


2.1 Two-level DP-matching algorithm 


The two-level DP-matching (TLDPM) algorithm’ attempts to find 
the best concatenation of reference patterns to match a given test 
pattern (of length M frames) by first determining the optimal reference 
pattern to match any portion of the test pattern and then attempting 
to find the optimal way in which to concatenate these pieces. Figure 
2a illustrates the first stage of this procedure. For all possible beginning 
points, 6, and all possible references, an isolated word dynamic time 
warping (DTW) algorithm is used to find the best path to all of the 
possible ending frames, F,. Distance scores for the best paths and the 
reference patterns which generates these best paths are recorded. The 
region over which a partial isolated word dynamic time warp is 
examined is shown in Fig. 2b. Here we show a region where the slope 
of the time warping function is restricted to be between 4 and 2 and 
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Fig. 2—Illustration of the TLDPM algorithm. 


where the maximum amount that the time alignment contour may 
deviate from the line of slope one, which starts at the point (0, 1), is 
plus or minus R frames. 

If we denote the accumulated distance for the uth reference pattern 
from starting frame } to ending frame e as Div, b, e), then the output 
of the first stage is the array of D values for all v, b, and e combinations. 
For any set of values of 6 and e, we can solve for the best reference 
pattern and distance as 


D(b, e) = min [D(v, 8, e)], (1a) 


N(b, e) = argmin [D(y, 8, e)], (1b) 


where D(b, e) is the minimum accumulated distance from frames 6 to 
e of the test pattern, and N(b, e) is the index of the reference pattern 
giving the minimum distance. 

The second stage of the TLDPM algorithm is to determine the best 
match by piecing together the reference patterns in an optimal manner. 
This is accomplished again using a dynamic programming algorithm 
which finds the best concatenation of / reference patterns, ending at 
frame e of the test pattern by trying all concatenations of a portion of 
the test pattern ending at frame e and all best reference pattern 
concatenations of length / — 1, i.e. 
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Di(e) = min [D(s, e) + Du-y(b — 1)] (2a) 


Do(0) = 0, (2b) 
where D,(e) is the accumulated distance generated by matching the 
best concatenation of / reference patterns to the portion of the test 
pattern between frame 1 and frame e and where D(0, e) is the distance 
associated with the best reference pattern used to match the portion 
of the test pattern between frame b and frame e. After computation of 
D,(e), the best string is recovered by first finding that value of J which 
minimizes D,;(M) and then tracing back through the sequence of 
decisions which were used to generate D;(M). 

The TLDPM algorithm made no inherent attempt to modify its 
reference patterns to account for either coarticulation or word shorten- 
ing. While word shortening is, in general, compensated directly by the 
DTW, the use of reference patterns which are inherently longer than 
the test patterns to which they are to be matched makes this problem 
much more difficult. In addition, the lack of any boundary modifica- 
tions for the reference patterns must inherently reduce the potential 
accuracy of the TLDPM approach. 

It is obvious that the TLDPM approach is not a sequential decision 
method, since no partial match is firmly decided on until the entire 
second pass of the algorithm is completed. It should also be clear that 
it is possible to generate not only the best string but the K best strings 
by simply keeping track of the K best reference patterns which match 
any portion of the test pattern and also by keeping track of the K best 
strings for every step of the second pass. Finally, it is clear that, as 
described, the TLDPM algorithm cannot make use of syntactical con- 
straints directly but must make use of them in a post-processing stage 
to choose among various candidate strings. That is, syntactical con- 
straints cannot be used to restrict which words of the vocabulary are 
used for any given beginning and ending point pair, unlike the LB 
algorithm in which levels correspond to word positions in the string. 
However, once the beginning and ending point pair distance matrices 
have been generated, then syntactical constraints can be applied in 
the second stage of the TLDPM algorithm. 


2.2 The sampling approach 


Unlike the TLDPM algorithm, the sampling algorithm is more of a 
sequential decision process. The sampling approach attempts, via a 
local minimum DTw algorithm,’’ to match a reference pattern to a 
portion of the test pattern. Unlike the two-level approach of the 
preceding section, not all portions of the test pattern are tested. 
Instead, only a small subset of the test pattern is used. The way in 
which the regions of the test pattern are chosen is as follows. Following 
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the time registration of the reference pattern to one portion of the test 
pattern, an ending point is implicitly defined by that frame of the test 
pattern which best matches the end of the reference pattern. After all 
reference patterns have been tried, the one which gives the best match 
is chosen as the proper word and the ending point of the test pattern, 
associated with this reference, is used to hypothesize a beginning 
region within the test pattern for the next set of references. This 
procedure is illustrated in Fig. 3. Reference pattern number 1, which 
is hypothesized to begin somewhere within beginning region 1, is time- 
registered to the test pattern using a local minimum DTW algorithm. 
Once the ending point for the best match for Ref. 1 has been deter- 
mined, a beginning region for the next reference pattern is determined. 
In general, this beginning region is centered somewhat earlier in the 
test pattern than the ending position of the previous word. This is 
done to compensate for the difference in durations between concate- 
nated reference patterns and connected word utterances, and to ac- 
count for coarticulation between words. 

Note that a categorical decision as to the best reference need not be 
made if the distance scores indicate a high likelihood of confusion. In 
such cases, when the distance scores for two or more reference patterns 
are approximately equal, the sampling method keeps tract of all 
possible strings using the approach described above. Thus, the possi- 
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Fig. 3—Illustration of the sampling algorithm. 
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Fig. 4—Two possible implementations of a constrained DTw algorithm. 


bility exists of finding a best string at the end of the search, which is 
different from the best string as the search proceeds from left to right. 

In addition to effectively modifying reference patterns by overlap- 
ping beginning and ending regions, the sampling algorithm can also 
use the local minimum DTW algorithm in such a way as to eliminate 
part of the end of any reference pattern. This is done by not forcing 
the local minimum Dtw algorithm to proceed all the way to the end of 
the reference pattern, but rather to stop at some frame before the end. 
Thus, the sampling approach has more inherent flexibility in dealing 
with modifications to the reference patterns than the TLDPM method. 

In contrast to the TLDPM algorithm, the sampling algorithm is able 
to use syntactical constraints (in the form of a regular grammar) 
directly, rather than in a post-processing stage. Such a method is 
described in detail by Levinson and Rosenberg,” but the main idea is 
to trace the graph which represents the grammar simultaneously with 
matching the reference patterns within the test pattern by keeping 
track of not only the current state, but also the current ending frame 
within the test pattern. 


2.3 The level-building algorithm 


Figure 4 illustrates the basic idea involved in the LB algorithm. Here 
we show how a constrained endpoint DTw algorithm, in which the 
slope of the warping function is restricted to be between % and 2, is 
used to find the best alignment between a test pattern and a given 
concatenated sequence of L reference patterns. In Fig. 4a, we show the 
computation proceeding in a sequence of vertical strips. Figure 4b 
shows an alternative way in which the computation may be performed. 
A set of horizontal lines has been drawn to indicate the boundaries 
between the different reference patterns in the concatenated reference 
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pattern. We may now compute the optimal, time alignment path in 
successive levels by using the distances accumulated at the end of one 
level to initialize the next level. The LB algorithm uses this decompo- 
sition to find the optimal sequence of reference patterns by trying all 
possible reference patterns at any level and recording, for each ending 
frame at that level, the reference pattern which gave the best distance 
to that ending frame and a pointer back to the previous levels. These 
minimum distances are used to initialize the following level. After the 
final level has been examined, the optimal path is recovered by tracing 
back along the chain of pointers. 

We have demonstrated that the algorithm solves exactly the same 
problem as the TLDPM algorithm.” In addition, we have shown how to 
modify the LB algorithm via a set of parameters to give it more 
flexibility. These parameters, as shown in Fig. 5, are as follows: 

(i) Sr, —Region of uncertainty at the beginning of the reference 
pattern. 

(ii) 45r, —Region of uncertainty at the end of the reference pattern. 
(tii) zeny —Region of uncertainty at the end of the test pattern. 
(tv) Mr —Multiplier used to reduce the size of the beginning region 

for any level. 

(v) | ¢€—Parameter used to restrict the size of any vertical strip. 
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Fig. 5—Ilustrations of the parameters of the LB algorithm. 
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Fig. 6—Generation of multiple candidate strings for the LB algorithm. 


The effects of these parameters are shown in Fig. 5. We observe that 
dr, and dz, define regions, at the beginning and end of each reference 
pattern, in which the local path may begin or end. In this manner, 
some of the gross features of coarticulation and length reduction 
present in a connected-word utterance may be accounted for. Thus, 
the LB algorithm has some of the flexibility inherent in the sampling 
algorithm. 

In addition to modifying templates, the algorithm has the advantage 
of not being a sequential decision process and, thus, has the ability to 
recover from mistakes as in the TLDPM algorithm. 

One important shortcoming of the LB algorithm is its ability to 
generate alternative candidates. The method by which multiple can- 
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didates are generated is to record not only the best candidate to each 
ending frame of each level, but to record the K best candidates and 
then to allow substitutions in the traceback. This method is illustrated 
in Fig. 6a for K = 2 candidates and for L = 2 levels. The solid paths 
represent the best paths from the beginning of the level to that 
particular ending point and the dashed paths represent the second 
best path to that particular ending point (using a different reference 
than the one used in the best path). We see that for K = 2 and L = 2 
we may generate four different candidate strings and, in general, may 
have K” different candidate strings. Such a situation may be repre- 
sented by a tree with a branching factor of K and a depth of L as 
shown in Fig. 6b for K = 3, L = 2. Generation of the possible candidate 
strings is simply a tree-searching problem, and it is possible to reduce 
the amount of traceback by pruning the tree. To prune the tree we 
may take advantage of the fact that the kth best path at any node 
in the tree represents a string whose distance is always larger than 
the (k — 1)st best path to that node. The difficulty with such a scheme 
is that, unlike the TLDPM algorithm, we are not guaranteed of finding 
even the true second best path. Figure 7 illustrates this problem. Here 
the best path is given by string AA and an alternative is given by string 
CA. Another path, shown by string BA, must have a larger distance 
than string AA, but may have a smaller distance than string CA. 
However, the path BA will not be recorded because, at the second 
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Fig. 7—Illustration of the failure of the LB algorithm to find all good candidate strings. 
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Fig. 8—Illustration of the path differences between the LB algorithm and a single time 
warp to the concatenated sequence of reference patterns. 


level, only the best path using reference A is recorded. While it would 
be possible to record this second best path, the computational problems 
which would have to be overcome seem much too burdensome. 

Another small theoretical difference between both the TLDPM and 
the LB solutions, and the exact best solution for the sequence of 
concatenated reference patterns that best matches the test pattern, is 
illustrated in Fig. 8. Here we show a sequence of two reference patterns, 
R, and Re, and the best path for warping the super reference pattern 
formed by concatenating R, and R, to T, as well as the LB path. Since 
the LB path is constrained to end at the end of each reference pattern 
(level), then a path point is constrained to occur at the end of each 
reference pattern. This need not be the case for the optimum con- 
catenation and, thus, a small difference can exist in the solutions. In 
practice, this effect does not occur (unless dr, = 5z,=0) since the 
5r,, Oz, parameters allow paths to skip frames of the reference patterns. 

A final consideration involving the LB algorithm is the use of 
syntactical constraints. Since all good candidates may not be generated 
by the LB algorithm, it is important that syntactical constraints be 
easily incorporated into its structure. Fortunately, this is possible and 
has been described in detail by Myers and Levinson.”* The basic 
principle is to build up the results by states of a finite-state automata, 
rather than by levels, and to use transitions of the finite state automata 
to guide the recognition process. 

In Section III, we summarize the qualitative features of the different 
connected-word recognition algorithms and also give a quantitative 
evaluation of both their space and time complexities. 
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lll. THEORETICAL COMPARISON OF THE CONNECTED-WORD 
RECOGNITION ALGORITHMS 


In Table I we summarize our qualitative evaluation of the different 
connected-word recognition algorithms. Both the sampling and LB 
algorithms have the ability to modify their reference patterns, although 
in different ways. The sampling approach allows both the removal of 

‘some frames from the reference pattern and the overlap of reference 
patterns, while the LB approach simply allows the removal of some 
frames from the reference patterns. Thus, we may conclude that the 
LB and the sampling algorithms will be able to more easily match test 
patterns which contain large amounts of coarticulation or length 
shortening than the TLDPM algorithm. 

We have shown that the sampling algorithm is a sequential decision 
process, while both the TLDpM and the LB algorithms have some form 
of backtracking. Thus, we expect better performance from the TLDPM 
and the LB algorithms in those cases in which the sampling procedure 
may get lost. 

In addition to being able to avoid getting lost, the TLDPM algorithm 
is the only algorithm that is capable of generating alternative candi- 
dates which are exactly correct. Thus, this algorithm should perform 
the best in cases in which many potential candidates are desirable. 
Unfortunately, as the final entry shows, the TLDPM algorithm is not 
well suited for syntactical analyses; therefore, it is more difficult to 
implement constraints using this algorithm.* 


Table |I—Qualitative comparison of connected-word recognition 


algorithms 
Modifi- 
cation of Sequen- Syntax 
Refer- _ tial Deci- Multiple Con- 
Algorithm ences sion Candidates straints 
Two-level DP-matching No No Exact No 
Sampling Yes Yes Heuristic Yes 
Level building Yes No Heuristic Yes 


3.1 Computational comparison 


In comparing both the time and space complexity of the three 
connected-word recognition algorithms, we shall assume that we are 
given a test utterance of length M frames, a vocabulary containing V 
words, a set of reference patterns with average length N (frames), with 
maximum time alignment deviation of +R frames. For purpose of 


* Note that it ts possible to implement both syntactical constraints and reference 
pattern modifications directly in the TLDPM algorithm without changing the fundamental 
ideas contained in the algorithm. For purposes of our comparison, however, we have 
considered the TLDPM algorithm only as it was originally described. 
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comparison, our measure of time will be the number of times that a 
frame of the test utterance is compared to a frame of any reference 
pattern. For the TLDPM algorithm, the number of comparisons is given 
by 


NCripem = V-M-N-(2R + 1), (3) 


since there is an isolated word time-warp of size N-(2R + 1) for all V © 
possible words, and for each of the M test frames. 
For the sampling algorithm, the number of comparisons is given by 


NC, = V-L-N.-(2€ + 1)-¥, (4) 


where L is the actual number of words in the test utterance, € is the 
range for the local minimum DTW algorithm, and jy is the average 
number of candidate strings which are retained. (Rabiner and Schmidt 
found that 7 was, on average, between 1 and 2)."’ Data in this figure 
result from one local minimum time-warp of size N -(2é + 1) for each 
possible reference, for each word of the test pattern and for each 
candidate string. 

For the LB algorithm, in which the slope of the warping function is 
restricted to be between % and 2 (as shown in Fig. 4), the number of 
comparisons is given by 


NCizs = V-Luax:-M-N/3, (5) 


because the size of parallelogram in Fig. 4 is about M-N-Lmax/3, and 
because all V references must be used at each level. 

Finally, by incorporating the range reduction techniques described 
by Myers and Rabiner,”* the number of comparisons required for the 
reduced LB technique becomes 


NCisr = V-Luax:-N-(2€ + 1), (6) 


since the reduced LB technique uses only a single local minimum time 
warp per level. 

Table IIa summarizes the computational aspects of these connected- 
word DTW algorithms. The row-labelled number of basic time warps 
refers to the number of times that a basic pTw algorithm is applied. 
The size of the time warp is the average size of these basic time warps, 
and the total is the number of comparisons given in eqs. (3) to (6). 

Table Ib gives a numerical comparison of the computation required 
by the various connected-word recognition algorithms for the case of 


Lomax = 5, V = 10, M = 120, N = 35, 
€=8c=12,R=12,7=15,L=4. 
Note that the total computation required by the TLDPM algorithm is 
15 times that of the LB algorithm and 30 times that of either the 
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Table lla—Computational comparisons of connected-word DTW algorithm 


Reduced Level 
Level Building Two-Level DP Matching Sampling Building 
Number of basic time warps Lax: V M-V L-V-y L.V 
Size of time warps N-M/3 _ N.-(2R + 1) N.-(2€ + 1) N.(2e€ + 1) 
Total computation for distance Imax: V+N-M/3 M.-.V-N-(2R + 1) L-V-y-N(2é + 1) L-V-N(2e + 1) 
Storage 3-M-Lmax:K 2-M.(2R+1)-K 0 3-M-Imax-K 


Table Ilb—Typical computational requirements for the case Luax = 
5,V=10,M=120,N = 35,€=8,€=12,y=1.5,L=4,K = 2, 


R=12 

Two-Level Reduced 
Level DP Level 

Building Matching Sampling Building 
Number of basic time warps 50 1200 60 40 
Size of time warps 1400 875 595 875 
Total computation for distances 70,000 1,050,000 35,700 35,000 
Storage 3600 12000 0 3600 


reduced LB or the sampling algorithm. The efficiency of the LB algo- 
rithm derives partly from the fact that optimal paths to a range of 
ending frames for a range of starting frames are found simultaneously. 


3.2 Storage comparison 


In comparing the storage required by the various connected-word 
recognition algorithms, we refer only to that storage which varies 
among the different algorithms. Hence, storage for reference patterns, 
and storage needed by all basic DTw algorithms is not considered here. 
The storage for the TLDPM algorithm is given by 


Stippm = 2-M-(2R + 1)-K. (7) 


Storage is required for the matrices of the K best distances and their 
associated words, for each of the M-(2R + 1) pairs of beginning and 
ending frames. 

The storage for the sampling algorithm is so small (less than 100) as 
to be negligible. The storage for the LB algorithm is given by 


Sip = 3-M-DImax-K (8) 


since, at the end of each level, it is necessary to store a distance, an 
associated word, and a pointer to the previous level for each of the K 
possible candidates. (Clearly, the factor M in eq. (8) can be substan- 
tially reduced if storage of distances, words, and back pointers is used 
only for finite distance ending frames.) 

The storage requirements for the four algorithms are summarized in 
the last row of Table Ila, and a numerical example (K = 2) is given in 
Table IIb. Note that both the LB and the TLDpM algorithm require a 
significant amount of storage, but that the storage required by the LB 
algorithm is only a third of that required by the TLDPM algorithm. 

We have seen that both the LB and the sampling algorithm are 
computationally simpler than the TLDpM, but that the TLDPM algo- 
rithm has the potential to generate better candidates for a given 
connected-word recognition task. Thus, certain trade-offs exist and 
some important questions must now be considered. Among them are 
the following: 

(t) Is the ability to modify the reference patterns a useful feature 
for a connected-word recognition algorithm? 
(zz) Is the ability to backtrack a useful feature, or is a sequential 
decision process sufficient? 

(uz) Are the additional computational costs of the TLDPM worth the 
increased number of good candidate strings? 

In Section IV we give the results of several experimental studies 
designed to answer these questions. 
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IV. EXPERIMENTAL COMPARISONS OF CONNECTED-WORD 
RECOGNITION ALGORITHMS 


To answer the questions at the end of Section III, each of the three 
connected-word recognition algorithms was simulated and tested on a 
common data base of 480 connected digit strings. The data base 
consisted of 80 strings of variable length (from two to three connected 
digits) spoken by each of six speakers (three male, three female). The 
strings were carefully articulated, and spoken in a slow, deliberate 
manner. Care was taken to guarantee an equal number of occurrences 
of all length strings, (from two to five digits), and an equal number of 
occurrences of all digits within the strings. This data base was used 
previously to evaluate the sampling method and the LB approach.’»”” 

Since results had been obtained previously on the sampling and LB 
algorithms, the TLDPM algorithm was all that had to be simulated and 
tested. For consistency, the Lpc feature set used in the sampling and 
LB systems was used in the TLDPM system as the basic analysis features. 
The only variable parameter in the basic two-level DP warp method is 
the time adjustment parameter R. Sakoe indicated that a value of R 
= 0.4 N would be appropriate where N is the average reference length.”° 
Since N = 40, this indicated that a value of R = 16 was required. The 
first experiment varied R from 8 to 99 to see its effect on recognition 
accuracy for a speaker-trained system. The results of this experiment 
(i.e., string error rate versus R) are given in Fig. 9. The results shown 
indicate that the larger the value of R, the lower the error rate, and 
that a flattening on the error curve occurs for values of R = 20. 


ERROR RATE IN PERCENT 
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Fig. 9—String error rate as a function of the time shift parameter R of the TLDPM 
method for connecting digit strings. 
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Fig. 10—-Percentage correct strings as a function of candidate position for the TLDPM 
method and the LB method for (a) no reference modification and for (b) reference 
modification. 


The next experiment made a direct comparison between the TLDPM 
method and the LB algorithm by setting the flexible LB parameters to 
values appropriate for the full range (« = 99, Mr = 99.9, Senp = 0, 
Tmiw = 9.9, Tmax = 9.9, Twin and Tmax specify distance thresholds 
used to reduce computation), and for no reduction in template lengths 
dr, = dz, = 0).”° The K = 2 best candidates were recorded at each level 
and the 10 best strings (of any length) were found from the LB output. 
Similarly, for the two-level DP warp, the 10 best strings were found 
directly. Figure 10a shows the percentage of strings correct as a 
function of candidate position in the list for both the TLDPM method 
(solid curve) and the LB method (dashed curve). Again, a speaker- 
trained system was used. The results show the same string error rate 
(11.5 percent) on the top candidate for both methods as anticipated 
from the earlier discussion. However, the results also show a 1 to 2 
percent higher string accuracy for the TLDPM method for candidate 
positions 2 to 10, indicating the potential improvements obtained by 
keeping distances for all posible beginning and ending frames. 
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Since the overall performance of the LB method was significantly 
better, using the best values of the LB parameters, rather than those 
shown in Fig. 10a, another experiment was run to compare recognition 
accuracy versus candidate position for both methods at the speaker- 
trained operating point (R = 99, e = 99, Mr = 1.4, denn = 4, Tun = 3.0, 
Tmax = 0.7, 5r, = 4, dr, = 6). It should be noted that in this case, the 
TLDPM algorithm was modified from its original specifications (to a 
form similar to the UE2-1 algorithm”’) to incorporate the parameters 
of the LB algorithm. These results are plotted in Fig. 10b. For the best 
candidate, the string error rates were 4.6 percent for the TLDPM 
method, and 4.8 percent for the LB method. For candidate positions 2- 
10, the TLDPM method has from 1 to 1.7 percent higher accuracy than 
the LB method. These results again indicate the small, but consistent 
improvement obtained by the TLDPM method, at least for alternative 
best string candidates. Whether or not this improved accuracy on 
higher candidate positions can be used depends strongly on the task, 
and the types of errors that occurred. Clearly, errors in string length 
can often be simply corrected, whereas errors of the same string length 
are generally difficult to detect and correct. 

A summary of the resulting string error rates on the top recognition 
candidate for all three connected-word recognizers is given in Table 
III. For the sampling and LB methods, results are also given for 
performance in a speaker-independent mode (using 12 isolated 
speaker-independent templates per digit). No results are given for the 
TLDPM system as a speaker-independent recognizer because of the 
inordinate amount of running time (about 4 to 5 hours/string) required 
to evaluate its performance. 


4.1 Performance of the LB method for other applications 


The LB method of connected-word recognition was also applied to 
the problem of recognizing names from a given directory from spoken, 
spelled input. Since the vocabulary here is letters of the alphabet, and 
because of the high degree of confusability among several of the letters 
(e.g., B, D, P, V, etc.), a somewhat different structure was used for 
recognizing the names.” First, a name class is found, and then all names 
(if any) within the name class are tested and a name score is stored. 


Table III—String error rates on connected digits 
Word Recognition Method 

TLDPM TLDPM 

(original) (modified) 


Speaker-trained 11.5% 4.6% 6.7% 4.8% 
Speaker-independent - - 9.0% 4.6% 


LB 
Sampling (reduced) 
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Table IV—Percentage names correct using the 


LB method 
Talking Speed 
Deliberate Normal 
Speaker-trained 99% 90.5% 
Speaker-independent 96% 87.5% 


The name with the smallest name score is chosen as the recognized 
name. 

To evaluate the performance of this system, four speakers (two 
male, two female) spoke 50 randomly selected names each as a con- 
nected sequence of letters of the last name, followed by a pause, 
followed by the initials. Each speaker spoke the name list two times. 
The first time, they spoke the names in a deliberate, carefully articu- 
lated manner, the second time, in a normal manner. The overall name 
accuracy achieved using the LB system is given in Table IV for both a 
speaker-trained and a speaker-independent implementation. Name 
accuracies of 96 to 99 percent were obtained for the deliberately spoken 
names, and of 87.5 to 90.5 percent for normally spoken names. 

In addition to directory assistance, the LB algorithm has also been 
applied to the problem of sentence recognition for an airlines reserva- 
tion system.’® In this experiment, four speakers (two male and two 
female) spoke 50 sentences each. These sentences were formed from a 
127-word vocabulary using a regular grammar. The grammar consisted 
of 144 states and 450 transitions. The sentences were spoken in a 
normal manner. Recognition was performed in a speaker-dependent 
manner. Sentence accuracies of 89 percent and word accuracies of 94 
percent were obtained. These experiments on directory assistance and 
sentence recognition demonstrate that the LB algorithm has wide- 
spread applicability. 


V. DISCUSSION AND SUMMARY 


The answers to two of the three questions at the end of Section III 
are clear. The gain in accuracy obtained by modifying the reference 
pattern is quite large for the connected-digit sequences (e.g., the string 
error rate falls from 11.5 percent to 4.6 percent for the TLDPM method). 
Similar improvements have previously been noted for the LB algo- 
rithm.’? Thus, it is clear that the modified reference patterns lead to 
improved overall performance. However, the results on recognizing 
spelled names (see Table IV) indicate clearly that as the rate of talking 
goes up, the performance of the system becomes dramatically worse. 
At this point, the basic model assumptions begin to fall apart and no 
simple modification of the reference patterns is adequate. 
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The answer to the second question concerning the value of a back- 
tracking versus a sequential approach can be obtained from Table III. 
By comparing performance of the sampling approach (with sequential 
decisions) to the LB or TLDPM results (with backtracking), we see a 
loss in string accuracy of 2 to 4.4 percent using the sequential approach. 
Since there is essentially no computational gain which accompanies 
this loss in accuracy, it seems clear that a backtracking approach is 
superior to a sequential decision method. 

The final question in section [V—whether the improved accuracy in 
higher candidate positions of the TLDPM method justifies the greatly 
increased computational load—is hard to answer. Since the recognition 
accuracy on the best string is the same, and since a cost factor in 
computation of about 40-to-1 is incurred, it would appear that the 
increased accuracy is not sufficient to justify the increased computa- 
tion. However, there may exist tasks with constrained syntax such that 
the improved accuracy does justify the costs. 

The overall conclusion is that both the LB and TLDPM methods 
provide high accuracy in recognizing strings of connected words. More- 
over, since the computation of the LB method is considerably less than 
the computation of the TLDPM, this method is to be preferred for most 
applications. 
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A Comparison of Three Speech Coders to be 
Implemented on the Digital Signal Processor 
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The recently developed digital signal processor is a device used for 
implementing low- to medium-complexity speech coders. It is cur- 
rently being used in implementing adaptive differential pulse-code 
modulation (ADPCM) coding, two-band sub-band coding, and four- 
band sub-band coding. This study was performed to determine opti- 
mal parameter values for the two sub-band coders in preparation for 
their implementation on the digital signal processor and to determine 
their performance relative to AappcM. (The actual implementation of 
the ADPCM and two-band sub-band algorithms are discussed in other 
papers in Part 2 of this issue of the Bell System Technical Journal.) 
Performance was judged on the basis of segmental signal-to-noise 
ratio and a forced-choice, subjective comparison test of the coders. 
All three coders were simulated at bit rates of 16, 20, 24, 28, and 32 
kb/s. The simulations were performed on a laboratory computer. 


I. INTRODUCTION 


The recently developed DsP is a device for implementing low-to 
medium-complexity speech coders. Three coders are currently being 
implemented. The simplest coder is adaptive differential pulse-code 
modulation (ADPCM) and is discussed in Ref. 1. The other two are in 
the sub-band coder (sBc) family. Of these, the simpler one is two-band 
sub-band coding (2B-sBc), featuring quadrature mirror filtering and 
two equal bands. It is discussed in Ref. 2. The other coder—the most 
complicated—is four-band sub-band coding (4B-sBc), featuring four 
equal bands. Its implementation is still in progress. 

This report discusses the initial design parameters for the latter two 
coders and the relative performance of all three. Segmental signal-to- 
noise ratio (SNR) measurements were made on all three coders via 
computer simulation at five different bit rates, 16, 20, 24, 28, and 32 
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kb/s. In addition, 12 subjects ranked the coders in a comparison test. 
The simulations reported here were carried out on a laboratory com- 
puter as preparation for the implementation of the sBc coders on the 
DSP. 

Section II reviews the design of ADPCM and discusses the design of 
the two sub-band coders. Section III discusses the results of the 
subjective testing experiment, and Section IV gives the conclusions of 
this study. 


ll. DESIGN OF THE CODERS 
2.1 Design of ADPCM 


The ADPCM design simulated here is based on the design of Cum- 
miskey et al.® A block diagram of the aDPcM coder described below is 
shown in Fig. 1. The most significant change from the design in Ref. 
3, is that only two multiplier values are used in changing the step-size, 
regardless of the bit rate. This is based on the ADPCM implemented by 
Johnston and Goodman.‘ This version of ADPCM has already been 
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Fig. 1—Adaptive differential pulse-code-modulation coder used in simulation. 
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implemented on the DspP and is described in Ref. 1. The ADPcM design 
was also used for quantizing the sub-band signals in the other coders. 

To simulate 20 kb/s and 28 kb/s ADPcM, alternating quantizers were 
used. For 20 kb/s, the two quantizers used are for 2 and 3 bits. Since 
the step-size is adapted based only on the most significant magnitude 
bit, the same step-size adaptation algorithm is used for all samples. 
The ratio of 2- to 3-bit quantizer step-size is held constant. This 
requires one additional multiplication to convert from the 2- to 3-bit 
step-size. 


2.2 Two-band SBC 


All sub-band coders are made from a few fundamental building 
blocks. The first is linear filtering to divide the signal into two or more 
sub-bands. These sub-bands can then be decimated to a lower sampling 
rate than the original signal. Some form of quantization must be used 
to encode and quantize each band. Interpolation and additional linear 
filtering is used to bring each band back to the original sampling rate 
and to its original space in the frequency spectrum. At this point, they 
can be added together to produce an output signal. 

The quadrature mirror filtering technique is fairly well known for its 
use with sub-band coders. Each pair of quadrature mirror filters (QMFs) 
produces two sub-bands of equal width in frequency. Johnson has 
compiled a collection of different length QMFs.” The possible quantizers 
which can be used are adaptive delta modulation (ADM), ADPCM, and 
adaptive pulse-code modulation (Apcm). Each of these techniques is 
fairly well known and understood. Likewise interpolation and deci- 
mation are also well understood. So what remains is the task of 
combining these building blocks in such a way as to fit on the DsP and, 
also, give the best possible performance. One of the tasks of this study 
was to choose good candidates for implementation. 

The 2B-sBC design is based on the 2B-SBc commentary grade coder 
of Johnston and Crochiere.® That coder was developed with the object 
of maintaining a high-quality Am radio signal. Its parameters were 
tuned to music rather than to speech. This section describes param- 
eters for a speech bandwidth version. There are five possible bit rates 
envisioned. For a more detailed discussion of sub-band coding in 
general and the exact implementation of this coder refer to Ref. 2. 

Figure 2 is a block diagram of the 2B-sBc. The input speech has been 
band-limited from 200 to 3200 Hz by a sharp bandpass filter and 
sampled at 8000 Hz. A 32-tap QmMF designed by Johnston is used for 
separating the digitized speech into the two sub-bands.” After 2-to-1 
decimation on both bands, we found average correlations for speech of 
0.7 and —0.45 for the low and high bands, respectively. 

The two bands are then coded using either ADPCM or ADM. The ADM 
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Fig. 2—Two-band sub-band coder used in simulation. 


is used only for the higher band at low bit rates. It is based on the ADM 
of Jayant.’ The prediction coefficients used for these coders were the 
correlation values mentioned above. Since the high band has a negative 
correlation, it was frequency inverted before quantization by the ADM, 
because this ADM requires a positive correlation for its adaptation 
mechanism to work properly. Since frequency inversion just means 
changing the sign of every other sample, this is a very minor operation. 

The next step was to determine optimal bit allocations for the low 
and high bands. After experimenting with different bit allocations and 
evaluating them on the basis of segmental SNR measurements and 
informal listening, the following bit allocations were adopted for the 
five bit rates: 


16 kb/s low: 3 bits high: 1 bit (ADM) 
20 kb/s low: 4 bits high: 1 bit (ADM) 
24 kb/s low: 5 bits high: 1 bit (ADM) 
28 kb/s low: 5 bits high: 2 bits 
32 kb/s low: 5 bits high: 3 bits. 


Some alternative designs were very close. For instance, a (4, 2) allo- 
cation for 24 kb/s is almost as good as (5, 1) for the speech it was 
tested on. Perhaps if the speech were less sharply bandpass-filtered 
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and if there were more high-frequency content (such as in telephone 
speech) the better allocation would be (4, 2) for 24 kb/s. 


2.3 Four-band SBC design 


The 4B-sBc design described here is new, although it is a logical 
extension of the 2B-SBc mentioned above. It starts with the same two 
sub-bands as the two-band design. Both of these bands are then 
divided into two new bands, yielding a total of four equally spaced 
bands. The filter used for the additional division in each band is the 
16-tap QMF of Johnston designated C in Ref. 5. Figure 3 shows a block 
diagram for this coder. 


TO 
MuLTi_ | CHANNEL 


PLEXER 





RECEIVER 


Fig. 3—Four-band sub-band coder used in simulation. 
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Once more the bands are quantized using ADPCM or ADM. Our 
measurements of average correlation for speech data showed correla- 
tions of 0.4, 0, 0, and 0.8 for the four bands going from lowest to highest 
in frequency. The fourth band (3000 to 4000 Hz) has actually been 
bandpass-filtered to cut off at 3200 Hz. As a result, it contains little 
power and can be ignored for low-bit rate coders. The correlations of 
the two middle bands are zero, reflecting that the long-term average of 
the speech spectrum from 1000 to 3000 Hz is flat. If a prediction 
coefficient of zero is used with ADPCM, the result is APcM. Thus, the 
two middle bands are AapcmM-encoded. The largest amount of power is 
in the first band; therefore, it receives the most bits. 

The bit allocations found to be the best by the same segmental SNR 
measurements and casual listening process were as follows: 


16 kb/s 4,2,2,0 (bands 1 to 4) 
20 kb/s 5,2,2,1 (ADM on band 4) 
24 kb/s 5,4,2,1 

28 kb/s 7,4,2,1 

32 kb/s 7,4,3,2. 


The greatest amount of error occurs in the lowest band. Even at the 
high rates (28 and 32 kb/s) this error is still perceptible as a low 
rumbling noise. However, it was found that a high-pass filter with a 
cutoff of 200 Hz eliminated this problem. The filter used was a 121-tap 
FIR filter. Table I gives the coefficients, and Fig. 4 shows the frequency 
response. A much smaller 11R filter could also be used to do the same 
job.” Note that the above bit assignments were made without using 
the FIR filter. With a high-pass filter, fewer bits could be allocated to 
the lowest band and more to bands two and three. 


2.4 Relative complexity of designs 


The appcm designs for 16, 24, and 32 kb/s have already been 
implemented on the psp.’ The combined encoder and decoder algo- 
rithms use 48 percent of the psp real-time capability for a sampling 


Table |—Coefficients for symmetric FIR high-pass filter 


—.153203E-01 —.488296E-03 —.427259E-03 —.305185E-03 —.152593E-03 
.610370K-04 .305185E-03 -610370E-03 .946074E-03 -131230E-02 
.170904E-02 .213630E-02 .259407E-02 .805185E-02 -350963E-02 
.3896740E-02 442518E-02 .485244E-02 .524918E-02 .558488E-02 
.589007E-02 -610370E-02 .625629E-02 .631733E-02 .628681 K-02 
.616474E-02 .592059E-02 .555437E-02 .509659E-02 448622E-02 

.387537BE-02 .283822E-02 -183111E-02 .732444E-03 —.549333E-03 

—.192267E-02 -—.344859E-02 —.506607E-02 —.677511E-02  —.857570E-02 

—.104679E-01 ~.123905E-01 —.144047E-01 —.164190E-01 —.184332E-01 

—.204779E-01 —.224921K-01 ~.244453E-01 —.263680E-01 —.281991E-01 

—.299387E-01 —.315867E-01 —.331126E-01 ~.344554E-01 —.356761E-01 

—.367443E-01 —.375988E-01 —.383007E-01 —.387890E-01 —.390942E-01 
.100000E-01 
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AMPLITUDE IN DECIBELS 
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FREQUENCY IN HERTZ 


Fig. 4—Frequency response of high-pass filter for 4B-SBC. 


rate of 8 kHz. An even lesser percentage of RAM and ROM memory is 
used. The 2B-sBc based on the design parameters reported here has 
also been implemented on the psp chip.” It uses 98 percent of the real- 
time capability and 78 percent of the RAM memory. It includes an 11R 
bandpass filter for the input. The 4B-sBc algorithm is planned for 
implementation in the near future. Since all of the major portions of 
the 4B-sBc have been programmed already for the 2B-SBc implemen- 
tation, it is possible to project how much of the DsP will be used. Both 
the transmitter and the receiver will require a DsP and each will use 
about the same fractions of real-time capability and RAM as the 
complete 2B-sBc algorithm. Therefore, so we might classify the three 
coders as having complexities of 0.5, 1, and 2, respectively. 


tl. RELATIVE PERFORMANCE OF THE CODERS 


Since the ssc designs are more complex, a demonstration of their 
improved performance over ADPCM was needed to justify their imple- 
mentation on the psp. To demonstrate their relative performance all 
three coders were simulated on a laboratory computer. Each processed 
speech from a stored file. The results were evaluated by both an 
objective and subjective measure. The objective measure was segmen- 
tal sNR, while the subjective measure was a forced-choice, subjective 
(A-B) comparison test in which all possible coders were compared: 

Six phonetically balanced sentences were used for evaluating the 
coders. Three were spoken by male speakers and three by females. 
They were recorded using a linear microphone. They were band- 
limited from 200 to 3200 Hz and sampled at 8000 Hz using a 15-bit 
linear quantizer. 


3.1 Segmental signal-to-noise ratio results 

In computing segmental SNR measurements, blocks of speech of 32 
ms were used. The ADPCM coder was compared with the original input 
speech. The sBc coders were compared with reassembled speech which 
had been processed by the appropriate QmF filtering, but with no 
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quantization. These slightly modified speech signals cannot be distin- 
guished from the original in casual listening. Without them it would be 
difficult to make a fair comparison of the three coders on the basis of 
SNR. The measurements on 4B-SBC were made before the 121-tap FIR 
high-pass filtering. 

The results of these measurements are summarized in Fig. 5. They 
show that the more complex SBC coders have a definite advantage over 
ADPCM at the lower-bit rates. Interestingly at 32 kb/s, ADPCM beats 
both of the more complex coders. The 4B-sBc maintains a fairly 
constant 2-dB advantage over 2B-SBC. In terms of bit rate this trans- 
lates to 4 kb/s. At the low rates, the 4B-sBc has about a 6-kb/s 
advantage over ADPCM. 


3.2 Subjective testing of the three coders 

An A-B comparison test was performed to rank the three coders. 
Each coder at each rate was compared twice against every other coder 
at every rate, as well as against the original. In the two comparisons of 
the two coders, each one was played in first position once. There were 
12 participants in the test and altogether there were 240 comparisons. 
The test was broken down into two parts, one with 110 comparisons, 
the other with 130. The participants listened over headphones in a 
soundproof booth. The participants were also broken down into two 
groups of six. If one group listened to a particular A-B comparison 
with a female speaker the other group heard a sentence with a male 
speaker and vice versa. Thus, we attempted to make a totally balanced 
and unbiased test. 


O ADPCM 
O 2B-SBC 
A 4B-SBC 


SIGNAL-TO-—NOISE RATIO IN DECIBELS 





16 20 24 28 32 
BIT RATE IN KILOBITS PER SECOND 


Fig. 5—Segmental SNR measurements for three coders. 
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Table Il—Coder versus coder ratings 


(Number represent percent obtained by coder listed at left of line against other coders) 


ADPCM 
Orig 16 20 24 
ADPCM 
16 0 = 17 0 
20 0 83 me 4 
24 0 100 9% — 
28 8 100 100 96 
32 2 100 100 #79 
2B-SBC 
16 0 67 21 8 
20 0 92 92 37 
24 4 100 96 83 
28 8 100 100° 87 
32 17 100 100 83 
4B-SBC 
16 8 100 92 25 
20 4 100 100 79 
24 4 100 96 79 
28 21 100 100 ~~ 92 
32 37 100 100 ~~ 87 


16 


33 
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96 
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63 33 
96 58 
42 13 
_— 29 
71 — 
87 54 
100 75 
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Fig. 6—Overall preference rankings for three coders. 


Table II gives the individual coder versus coder comparisons. In 
addition, an overall preference ranking was computed based on the 
total number of votes received by each coder. In all, a total of 360 
votes could be received by any coder. Figure 6 shows the percentage 
of the 360 possible votes received for each coder. This result is in good 
agreement with the results of Fig. 4. For example, both sub-band 
coders show adantage over ADPCM at the low rates and ADPCM catches 
up or passes them at the high rates. 

Some of the more significant results are the following: 

(t) The 4B-sBc has an 8-kb/s perceptual advantage over ADPCM at 
the low rates. The 24-kb/s appcM has been used for voice storage and 
playback systems.® This result indicates that 16-kb/s 4B-sBc could be 
substituted at a 33 percent savings in storage or, equivalently, a 50 
percent increase in message storage capability. Moreover, at 20-kb/s, 
2B-SBC has a 4-kb/s perceptual advantage over ADPCM. 

(ti) Although 4B-sBc lost to ADPCM at 32 kb/s in SNR measure- 
ments, it beat ADPCM in the subjective tests. In addition, in direct 
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comparisons with the original, 32-kb/s 4B-SBC received 37.5 percent of 
the votes, an almost equipreferential rating. This indicates it is high 
quality. Since 32-kb/s ADPcm is often described as toll quality, then 
32-kb/s 4B-SBC also deserves this ranking. 

(tit) The 2B-SBC seems to provide a good alternative to ADPCM at 
the low bit rates for a modest increase in complexity. 


IV. CONCLUSIONS 


We have presented measurement data and simulation results for use 
in implementing two sub-band coders on the pDspP. This data has already 
been used to design and implement the 2B-sBc and is being used for a 
planned 4B-sBc implementation. Simulations of these candidate coders 
were made on a laboratory computer. The results of these simulations 
indicate 2B-sBc and 4B-sBC have important advantages over ADPCM 
at low bit rates. This advantage is as much as 8 kb/s for 4B-sBc and 
4 kb/s for 2B-sBc. The 16-kb/s 4B-sBc could be substituted for 24-kb/ 
S ADPCM in a voice storage and playback system. In addition, 4B-sBc 
is rated as better quality at 32 kb/s than at 32-kb/s ADPCM. 

Since the complexity of these coders is within an order of magnitude 
of ADPcM they should be considered as viable alternatives. 
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